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SECTION  I 


Introduction 


Aircraft  flutter  is  a  destructive  phenomenon  which  requires  special  attention  in 
the  design  process.  The  elements  of  flutter  are  structural  dynamics  and  unsteady 
aerodynamics.  Of  these,  it  is  generally  recognized  that  unsteady  aerodynamics 
are  the  more  difficult  to  model  and  the  least  reliable.  In  1935,  Theodorsen  was 
the  first  to  develop  a  practical  unsteady  incompressible  aerodynamic  formula1 
for  a  flutter  analysis  of  a  two  dimensional  airfoil.  It  was  fifty  years  ago  that 
Smilg  and  Wasserman  of  the  Aircraft  Laboratory  of  the  Wright  Air  Development 
Center  wrote  their  landmark  report  on  flutter  clearance  using  the  K-method  and 
strip  theory.  Of  course  such  methods  can  be  addressed  with  manual  calculations. 
Compressibility  is  normally  associated  with  the  flutter  of  high  speed  aircraft  It 
is  impractical  to  solve  the  compressible  unsteady  aerodynamic  equations  by 
hand.  The  doublet  lattice  method2  was  developed  along  with  improvements  in 
digital  computer  technology.  Hopefully,  the  doublet  lattice  method  represents 
the  most  rudimentary  unsteady  aerodynamic  technique  in  practice  where 
subsonic  compressible  flow  is  a  consideration.  With  the  introduction  of  today’s 
supercomputers,  non-linear  aerodynamics  are  now  being  addressed,  in  spite  of 
the  high  cost.  It  is  because  of  the  high  cost  vxhnical  complications 
associated  with  non-linear  Computational  Fluid  Dynamics  (CFD)  that  the 
doublet  lattice  method  is  still  used  almost  exclusively  for  die  subsonic  flutter 
clearance  of  flight  vehicles  being  designed  today.  It  is  difficult  to  imagine  the 
day  when  non-linear  CFD  will  replace  the  doublet  lattice  method  in  the 
preliminary  design  environment 


1.  Section  $-6  of  the  wade  by  Bijplinjboff  m.  &L  ooot tint  «n  MpUnalioo  of  Tbeodoktto'i  fottmiU. 

2.  GiMtOf.  Kalaua  end  Rotten 
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While  this  document  is  not  a  survey  report,  it  is  appropriate  to  acknowledge  the 
original  authors  of  the  doublet  lattice  method,  Dr.  Edward  Albano  and  Dr. 
William  P.  Rodden.  The  subsequent  work  of  Mr.  Joseph  P.  Giesing,  Mrs.  Terez  P. 
Kalman,  and  Dr.  William  P.  Rodden  of  the  Douglas  Aircraft  Company  was 
sponsored  by  the  Air  Force  Flight  Dynamics  Laboratory  under  the  guidance  of 
Mr.  Walter  J.  Mykytow.  The  two  computer  codes  which  resulted  from  this 
contractual  effort  are  H7WC1  and  following  that,  N5KA2.  These  codes  are  still 
the  de  facto  standard  where  the  flutter  clearance  of  military  aircraft  is  involved. 
The  geometric  options  offered  in  these  codes  are  extensive,  including  multiple 
surfaces  and  slender  bodies.  The  purpose  of  this  document  is  to  derive  the 
fundamental  formulae  of  the  doublet  lattice  method.  In  order  to  keep  focused  on 
the  fundamentals,  the  formulae  derived  in  this  report  are  restricted  to  planar 
wings.  The  additional  work  to  extend  the  formulae  to  wings  with  dihedral  is  not 
conceptually  significant  Unsteady  aerodynamics  over  slender  bodies  is  not 
addressed  here. 

The  mathematical  background  leading  to  the  doublet  lattice  method  is  found 
among  many  documents  and  texts.  Considering  the  importance  of  the  doublet 
lattice  method,  it  seems  surprising  that  we  lack  a  single  consistent  derivation. 
This  document  attempts  to  answer  the  need  for  a  unified  derivation  of  all  the 
important  formulae  from  first  principles  to  the  integral  formula  and  also  includes 
a  simple  doublet  lattice  source  code.  The  target  audience  is  the  graduate  student 
or  engineer  who  has  had  a  first  course  in  aeroelasticity  and  would  like  to  focus 
on  the  mathematics  of  the  doublet  lattice  method.  The  author  has  assumed  the 
reader  has  a  familiarity  with  the  classical  topics  of  potential  aerodynamics  and 
linear  boundary  value  problems. 

The  author  takes  no  credit  for  developing  the  formulae.  All  the  work  presented 
here  was  compiled  from  many  references  to  create  a  unified  derivation.  The 
author  does  take  credit  for  any  additional  illumination  which  he  may  cast  on 


1 .  Gtesing  el  *1,  AFFDL-7 1-5.  Vo!  EL.  Pan  ]  it  the  origin*]  pilot  code,  ii  u*e*  doublet  panel*  to  model  bodie*  u  uutultr  wing*. 

2.  See  Gieiing  et  *1.  AfFDL-7 1-3.  Vo]  H  Put  U  i*  tbe  fin*]  deliverable  code  tod  uae*  *x>*i  doublet*  sad  interference  panel* 
for  bodie*  u  well  as  the  method  of  images. 
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these  derivations.  The  main  contribution  of  this  report  is  that  all  these 
derivations  are  presented  in  a  logical  sequence  in  a  single  document  not 
available  elsewhere.  The  hope  is  that  the  reader  will  gain  an  accurate 
appreciation  of  the  doublet  lattice  method  by  following  this  single  derivation. 

This  report  focuses  on  presenting  the  general  mathematical  procedure  behind  the 
doublet  lattice  method.  Such  mathematics  do  not  make  easy  reading.  The  task  of 
making  these  mathematical  derivations  pleasurable  may  be  impossible.  Learning 
these  mathematics  requires  that  one  take  pen  and  paper  in  hand  and  derive 
unfamiliar  formulae.  The  integral  formulae  of  Sections  XII  and  XIII  may  seem 
excessively  complicated.  However,  this  complication  is  a  matter  of  bookkeeping 
and  not  a  matter  of  high  level  mathematics  beyond  undergraduate  calculus. 

In  short,  the  doublet  lattice  method  is  based  on  the  integral  equation  (276).  The 
integrand  of  this  equation  models  the  effect  of  the  pressure  difference  (across  the 
plane  of  the  wing)  at  one  wing  location  on  the  induced  upwash  (component  of 
velocity  which  is  normal  to  the  plane  of  the  wing)  at  another  wing  location.  In  a 
sense,  it  can  be  said  that  equation  (276)  is  entirely  equivalent  to  the  linear 
aerodynamic  potential  equation  (42)  and  the  linearized  pressure  equation  (52). 
While  equation  (276)  is  a  specialization  of  equations  (42)  and  (52),  no 
Approximations  were  assuramed  in  its  derivation  from  the  Euler  equations. 

Euler’s  five  differential  equations  of  in  viscid  flow  are  the  starting  basis  for  all 
derivations  in  this  report  These  five  equations  are  comprised  of  one  equation  of 
continuity,  three  equations  of  momentum  and  one  equation  of  state.  The 
equations  of  momentum  model  pressure  equilibrium  in  each  of  three  coordinate 
directions.  The  inviscid  restriction  of  Euler’s  equations  means  the  momentum 
equations  lack  terms  of  shear  force.  With  no  shear  force  on  a  fluid,  no  vorticity 
(flow  rotation)  can  be  developed.  While  the  Euler  equations  are  restricted  from 
generating  rotation,  they  are  not  restricted  from  convecting  rotation  if  rotation 
exists  in  the  initial  or  boundary  conditions  to  Euler’s  differential  equations.  This 
is  the  starting  assumption  in  all  the  subsequent  mathematical  developments. 
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The  solution  to  a  boundary  value  problem  satisfies  both  the  coupled  partial 
differential  equations  and  the  associated  boundary  condition  equations.  When  a 
single  unique  solution  exists,  the  number  of  variables  equals  the  number  of 
partial  differential  equations.  For  the  Euler  equations,  we  have  five  variables  and 
five  differential  equations.  These  five  equations  describe  the  flow  within  the 
domain.  The  boundary  condition  is  a  description  of  the  flow  on  the  boundary  of 
the  domain.  The  important  point  here  is  that  once  Euler’s  boundary  value 
problem  has  been  stated,  all  that  remains  is  the  mathematical  solution.  After  a 
general  mathematical  procedure  has  been  identified,  engineers  can  proceed  to 
make  automated  applications  to  their  design  procedures. 

Section  13  starts  with  Euler’s  differential  equations  with  five  unknown  variables 
representing  flow  density,  pressure  and  three  components  of  velocity.  Euler’s 
equations  are  non-linear.  The  doublet  lattice  method  is  linear.  The  objective  of 
Section  II  is  to  reduce  Euler’s  non-linear  boundary  value  problem  (five  differen¬ 
tial  equations  and  associated  boundary  conditions  with  five  unknown  variables) 
to  a  linear  boundary  value  problem  (one  linear  differential  equation  (42)  with 
linear  boundary  conditions  in  terms  of  one  unknown  variable,  the  velocity 
potential).  The  derivation  of  this  linear  (potential)  equation  follows  the  approach 
taken  in  the  text  by  Bisplinghoff,  Ashley  and  Halfman.  The  text  by  Karamcheti 
provides  an  excellent  explanation  of  the  velocity  potential.  Equation  (42)  and  the 
boundary  conditions  developed  in  Section  IV  are  sufficient  for  generating  a 
unique  solution  for  the  velocity  potential. 

Section  III  assumes  the  velocity'  potential  of  Section  II  is  now  a  known  quantity. 
The  objective  of  Section  HI  is  to  develop  a  single  linear  differential  equation  for 
the  unknown  pressure  variable.  After  all,  the  aerodynamicist  is  interested  in  the 
pressure  loads  on  the  aircraft,  not  the  velocity  potential.  The  desired  linear  for¬ 
mula  is  equation  (52).  The  remainder  of  Section  III  provides  the  derivation  of 
die  reverse  relation.  In  other  words,  given  the  pressure  over  the  domain,  what  is 
the  potential  over  the  domain.  This  relation  is  equation  (73). 


4 


Introduction 


If  one  accepts  the  linear  differential  equations  (42)  and  (52)  as  a  linear  model  for 
the  flow,  then  Euler’s  non-linear  differential  equations  are  somewhat  irrelevant 
to  the  subsequent  mathematical  development  of  the  doublet  lattice  method. 
Again,  all  that  remains  is  the  identification  of  a  solution  which  satisfies  the  aero¬ 
dynamic  potential  equation  (42)  and  the  boundaiy  condition. 

In  Section  IV,  the  non-linear  tangential  flow  boundaiy  condition  for  Euler’s 
boundary  value  problem  are  stated  and  then  linearized  The  linear  boundaiy 
condition,  together  with  the  linear  aerodynamic  potential  equation  (42)  form  the 
boundary  value  problem  that  will  be  solved  by  the  doublet  lattice  method. 
Finally,  the  special  case  of  the  boundary  condition  on  an  oscillating  wing  is 
presented.  The  form  of  the  doublet  lattice  method  presented  here  assumes  the 
wing,  and  therefore  the  flow,  oscillate  harmonically.  Complex  notation  is 
assumed  here  and  the  reader  must  be  familiar  with  solving  complex  algebra 
problems. 

At  this  point,  all  the  preliminary  aspects  leading  to  the  doublet  lattice  method 
have  been  completed.  The  linear  boundary  value  problem  has  been  completely 
described.  The  solution  procedure  begins  to  take  shape  in  Section  V. 
Mathematicians  will  typically  solve  simple  linear  boundary  value  problems 
using  the  method  of  separation  of  variables.  This  approach  is  not  at  all  practical 
for  solving  the  flow  over  even  simple  wing  planforms.  Another  method  is  to 
identify  a  set  of  solutions  to  the  linear  aerodynamic  potential  equation  (42).  If 
these  solutions  can  be  linearly  superimposed  such  that  the  boundary  conditions 
are  at  least  approximately  satisfied,  then  the  solution  is  complete. 

Simple  solutions  to  the  aerodynamic  potential  equation  are  not  easily  identified. 
This  is  the  motivation  for  Section  V  in  which  we  transform  the  aerodynamic 
potential  equation  to  the  well  known  acoustic  potential  equation.  We  will 
identify  an  elementary  source  solution  to  this  acoustic  equation  in  Section  VI. 
Finally,  this  solution  is  modified  and  transformed  back  to  the  coordinates  of  the 
aerodynamic  potential  equation  in  Sections  VI,  VII  and  Vffi,  Equation  (152)  is 
the  elementary  source  solution  (x,  y,  s*  4,  T),  t)  to  the  aerodynamic 
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potential  equation.  The  definitions  for  R  and  x,  found  in  equation  (152)  are 
given  in  equations  (149)  and  (151).  The  arguments  x ,  y, :  and  t  in  the  function 
<j>  are  identical  to  the  arguments  of  the  potential  <f>  (x,  y,  z,  t ) .  The  arguments 
r\  and  C,  represent  the  a\  y  and  z  coordinates  of  the  reference  point  which  labels 
the  elementary  source  solution.  In  other  words,  each  source  solution 
<j>s  (x,  y,  z,  r\,  t)  represents  the  potential  at  coordinates  (x,  y,  z,  t )  due  to  a 
spherically  symmetric  flow  originating  from  a  point  located  at  4,  "H  and  C- 

Clones  of  the  elementary  source  solutions  can  be  placed  throughout  the  flow 
domain.  The  resulting  potential  field  is  a  superposition  of  the  potential  arising 
from  each  point  source  as  given  by  equation  (152).  It  turns  out  that  the  point 
sources  are  placed  at  coordinates  T),  £  on  the  boundary  of  the  vehicle.  If 
insufficient  point  sources  are  used,  the  approximating  composite  potential 
solution  will  be  “bumpy”.  In  order  to  provide  sufficient  smoothness,  the 
approximating  solution  requires  sufficient  point  sources  on  the  surface.  As  a 
matter  of  fact,  one  can  take  the  limiting  case  of  a  continuum  of  point  sources  on 
the  surface.  Each  differential  area  of  the  aerodynamic  boundary  (excluding  the 
far  field  boundary)  contains  a  coordinate  point  identifying  another  clone  of  the 
elementary  solution  with  its  own  strength.  The  linearized  aerodynamic  boundary 
on  a  wing  is  a  two  dimensional  plane  or  sheet. 

Section  IX  describes  a  continuous  source  sheet  A  source  sheet  is  a  two 
dimensional  surface  which  has  been  partitioned  into  differential  areas.  Each 
differential  area  is  assigned  a  point  sour'?  of  varying  flow  rate  (or  strength).  In 
the  limit  as  the  partition  is  infinitely  refined,  a  continuous  source  sheet  is 
formed.  This  purpose  of  this  section  is  for  illustration  only.  A  source  sheet  will 
not  solve  all  the  boundary  conditions  for  flow  over  a  wing.  The  potential  field 
associated  with  a  source  sheet  is  the  same  above  and  below  the  plane  of  the 
sheet  The  pressure  field  is  entirely  dependent  on  the  potential  field.  If  the 
potential  is  the  same  above  and  below,  then  it  is  not  possible  for  a  single  source 
sheet  to  generate  a  pressure  difference.  However,  two  opposing  source  sheets 
can  generate  a  pressure  difference. 
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Section  X  describes  a  continuous  doublet  sheet  This  is  the  limiting  condition  as 
two  source  sheets  are  brought  together.  Each  sheet  has  an  opposite  strength 
proportional  to  the  inverse  of  the  distance  between  them.  A  single  doublet  sheet 
can  be  used  to  mathematically  model  the  pressure  difference  between  the  upper 
and  lower  surfaces  of  a  thin  wing.  One  could  use  this  doublet  formulation  to 
solve  a  boundary  value  problem  for  flow  over  a  thin  wing.  However,  this  doublet 
sheet  formula  is  in  terms  of  the  (velocity)  potential.  We  are  really  interested  in 
the  pressure  load.  Therefore,  given  the  solution  for  the  potential,  we  then  have  to 
solve  a  second  problem  using  differential  equation  (52)  to  obtain  the  pressure 
field  from  the  potential  field.  A  more  direct  approach  is  taken  in  the  remaining 
sections.  Unfortunately,  this  direct  approach  increases  the  complexity  of  the 
formulation. 

In  Section  XI,  the  pressure  potential  and  the  acceleration  potential  are 
introduced.  The  pressure  potential  is  the  unknown  variable  of  the  pressure 
potential  equation  (185).  The  pressure  potential  equation  arises  as  a  direct 
consequence  of  the  aerodynamic  potential  equation  (42)  and  the  pressure 
formula  (52).  The  form  of  the  pressure  potential  equation  is  identical  to  the 
aerodynamic  potential  equation.  Therefore,  any  elementary  solution  to  the 
aerodynamic  potential  equation  is  also  a  solution  to  tire  pressure  potential 
equation.  Therefore,  the  elementary'  solution  to  equation  (185)  is  a  source 
function.  Again,  the  pressure  source  sheet  is  symmetric  and  cannot  generate  a 
pressure  difference  above  and  below'  its  plane.  A  pressure  doublet  solution  for 
oscillatory  flow',  otherwise  known  as  the  acceleration  potential  is  developed 
in  equation  (191).  A  pressure  doublet  sheet  can  generate  a  pressure  difference. 
Now  that  an  elementary'  pressure  doublet  solution  has  been  identified,  it  is 
necessary'  to  formulate  the  potential  field  (and  then  the  velocity  at  the  wing 
boundary)  which  arises  with  this  pressure.  This  is  given  in  equation  (197). 

In  Section  XII,  we  use  the  point  pressure  doublet  solution  of  Section  XI  to  create 
a  pressure  doublet  sheet  using  the  same  approach  taken  earlier  to  expand  the 
point  source  to  a  source  sheet  This  results  in  integral  equation  (203)  which 
describes  the  acceleration  potential  (pressure)  which  arises  from  a  pressure 
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doublet  sheet  What  we  really  need  and  subsequendy  derive  is  Equation  (224) 
which  is  the  w  component  of  velocity  which  arises  with  a  pressure  doublet  sheet 
This  is  the  integral  formula  of  Section  XII.  If  one  determines  a  pressure  doublet 
distribution  which  satisfies  the  tangential  flow  boundary  conditions  using  this 
integral  formula,  then  the  boundary  value  problem  has  been  solved.  In  a  sense, 
the  linear  boundary  value  problem  and  the  specialized  integral  formula  are 
equivalent  While  the  concept  is  simple  enough,  the  procedure  for  determining 
the  distributed  strength  of  the  pressure  doublet  is  not  immediately  obvious. 

Within  the  integral  formula  is  the  kernel  function.  The  kernel  function  is  highly 
singular  (division  by  zero  valued  variables  of  power  higher  than  one)  at  the 
surface  of  the  pressure  doublet  sheet.  The  purpose  of  Section  XHI  is  to  reduce 
the  severity  of  the  singularity  and  to  put  the  kernel  function  in  a  form  which 
lends  itself  to  numerical  evaluation  with  a  computer.  (Note:  the  singularity  of  the 
integral  formula  is  not  entirely  removed.)  Section  Xm  is  very  detailed  and 
inherently  difficult  to  follow  with  all  the  variable  substitutions.  The  kernel 
function  for  a  planar  wing  is  simply  summarized  in  the  following  section  as 
equation  (277)  with  supporting  equations  (278)  through  equation  (281). 

The  doublet  lattice  method  is  a  solution  procedure  for  the  integral  formula.  The 
normal  velocity  component  w  is  the  known  boundary  condition.  What  we  don’t 
know  is  the  pressure  difference  Ap  across  the  thin  wing.  Unfortunately,  the 
unknown  pressure  falls  within  the  integrand.  In  the  doublet  lattice  method,  the 
pressure  doublet  sheet  is  divided  into  trapezoidal  areas.  Within  each  trapezoid, 
the  unknown  pressure  function  is  assumed  to  take  a  form  with  unknown  constant 
coefficients.  More  specifically,  the  pressure  is  assumed  spatially  constant  within 
the  trapezoidal  area.  The  integral  formula  is  evaluated  over  each  trapezoidal  area 
independently.  The  unknown  constant  coefficients  come  out  from  under  the 
integrand.  Some  of  the  remaining  singularities  in  the  integrand  are  avoided  by 
replacing  the  pressure  doublet  sheet  with  a  pressure  doublet  line.  Remaining 
singularities  in  the  resulting  line  integul  (over  the  doublet  line)  are  addressed 
using  the  concept  of  principle  values. 
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The  formulae  of  Section  XIV  are  encoded  in  the  doublet  lattice  program  of 
Appendix  A  as  described  in  Section  XV.  The  reader  must  keep  in  mind  that  the 
solutions  obtained  here  will  not  agree  exactly  with  the  results  of  Giesing  et.  al. 
because  in  this  report,  there  have  been  no  steady  state  corrections.  Furthermore, 
the  formulae  of  Section  XIV  are  restricted  to  a  single  wing  in  the  (x,  y)  plane. 

One  misconception  about  the  doublet  lattice  method  is  that  the  formulae  are 
made  non-singular  by  replacing  the  doublet  sheet  with  a  lattice  structure.  This  is 
not  the  case.  The  doublet  lattice  method  is  made  non-singular  only  by  using  the 
concept  of  principle  values.  The  appeal  of  the  doublet  lattice  method  is  solely  in 
its  programming  simplicity.  The  assumed  form  of  the  pressure  function  would 
normally  change  according  to  the  proximity  to  the  edges  of  the  wing.  With  the 
doublet  lattice  method,  die  form  is  the  same  for  all  elements.  This  is  a  major  sim¬ 
plification. 

In  providing  this  overview  of  the  doublet  lattice  method,  the  reader  will 
hopefully  be  better  prepared  to  investigate  the  mathematical  details  of  the 
following  report.  Again,  the  interested  reader  will  find  that  by  using  pen  and 
paper  while  reading  this  report,  he  will  obtain  a  special  level  of  ownership  of  the 
material  presented  here.  There  leally  is  no  other  way  to  learn. 
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From  First  Principles  to  the  Linearized  Aerodynamic  Potential  Equation 


The  equations  describing  an  inviscid  fluid  flow  over  a  solid  body  will  be 
described  in  a  frame  of  reference  that  is  attached  to  the  body  and  travels  with  it 
Here,  the  three  cartesian  components  are  indicated  by  x,  y  and  z.  Time  t  is  a 
fourth  coordinate  or  dimension.  The  five  state  variables  are  pressure, 
p  (x,  y,  z,  t) ,  density,  p  (x,  y,  z,  t) ,  and  the  three  cartesian  components  of  the 
velocity  with  respect  to  the  frame,  u  (x,  y,  z,  t) ,  v  (x,  y,  z,  t) ,  and  w  (x,  y,  z,  t) . 
A  control  volume  is  identified  and  five  appropriate  equations  are  formulated  to 
solve  for  the  five  state  variables.  These  are  the  partial  differential  equation  of 
continuity,  three  partial  differential  equations  of  momentum  and  an  equation  of 
state.  These  formulae  are  presented  here  without  derivation1. 

The  continuity  equation  is  given  as 


dp  d(p«)  d(pv)  d(pw) 

___  +  — «- — .  — «- —  -t-  — «- —  —  u 

at  dx  dy  dz 


(1) 


The  three  components  of  the  momentum  equations  are  given  here.  For  momen¬ 
tum  in  the  x,y  and  z  directions 


du  du  du  9  u 

-1 3P 

(2) 

*sr—  4-21^—  = 

at  dx  ay  dz 

p  < h 

dv  dv  dv  dv 

-ldP 

(3) 

K-  4-ttvr-  +  Vvr-  -t-  W*~  = 

at  dx  dy  dz 

JTy 

1.  One  ii  directed  to  either  Ohjkct  5  of  Kuraorcbeti  o»  Appt.xiix  3  of  Koethe  inH  Chow,  lo  the  Utter,  th-s  full  Navier-Stdcee 
etjuetiaM  aue  derived  from  which  the  Euler  cqutiioai  ire  obaiaed  by  Mihnjt  the  'Ucouitereu  to  texo. 
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3w 

u 


3w 


+  w 


dw  _  -1  dP 
Tz  -  -pTz 


(4) 


With  the  following  isentropic  relation1,  we  complete  the  desired  set  of  five  equa¬ 
tions  which  we  use  to  solve  for  the  five  unknown  state  variables  u,  v  w,  p  and 

P- 


(5) 


The  constant  y  is  the  ratio  of  specific  heats.  The  variables  p0  and  po,  are 
constant  referer  t  values  of  pressure  and  density.  This  ratio  of  pressure  and 
density  is  constant  for  any  element  of  fluid.  If  the  entire  fluid  field  originates 
from  a  single  static  reservoir  then  the  ratio  of  equation  (5)  is  constant  for  the 
entire  fluid  field.  The  formula  for  the  apeed  of  sound  is  also  presented  here 
without  development2. 


a 


2  dp 
dp 


As  a  first  step  in  reducing  the  problem,  we  denne  die  fluid  velocity  vector 


(6) 


q  a  tfi4v)  +  4  ^ 

The  condition  of  irrotationality,  Vxq  »  0,  allows  us  to  introduce  the  velocity 
potential  <S>  (x,  y,  z,  k) ,  a  new  state  variable,  the  relationship  between  the  veloc¬ 
ity  potential  and  the  three  unknown  velocity  components  is  presented  here  with¬ 
out  further  explanation  \ 


1 .  Equation  (5)  it  the  cLimuJ  aqairioo  of  «Ufe  restricted  to  iaotrook  condition!.  Such  fluid  flow  U  called  bwotropk.  See 
page  50  of  Liepmaon  end  Rothko. 

2.  S'*  rhiptcr  L8  of  Lkpjui^f.  aud  Rothko.  Alternative]?,  aee  chapter  92  of  Kuetbe  and  Chow. 

5.  See  page  244  of  Kaiamcbeti. 
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VO  = 


.do  .do  .do 

,^+J'3 y+kH 


(8) 


With  the  velocity  potential,  we  can  reduce  our  five  equations  and  five  unknowns 
to  three  equations  and  three  unknowns.  The  three  unknowns  are  O,  p  and  p. 
These  three  equations  will  be  (23),  (18)  and  (5). 


We  denote  the  magnitude  of  the  velocity  as  q.  It  should  now  be  clear  from  equa¬ 
tions  (7)  and  (8)  that  the  relationship  between  the  velocity  magnitude  and  the 
velocity  potential  is 


q  =  Ju2  +  v2  +  w2  =  J<P2  +  $2  +  <X>2  (9) 

We  can  explicitly  state  the  basis  for  the  irrotational  condition  V  x  q  =  0  in 
terms  of  the  following  three  relations.  This  is  obtained  using  the  components  of 
q  given  in  equation  (7)  or  equation  (8). 


dw>  dv 

dy  ~  5z 

dw  du 
dx  ”  dz 


(10) 

(11) 


dv  du 
Tx  =  5 


(12) 


The  remainder  of  this  section  elaborates  on  chapter  5.1  of  the  work  by  Bispling- 
hoff  Ashley  and  Halfman  (or  Dowell  et.  al.).  The  three  momentum  equations  (2) 
-  (4)  can  be  put  into  a  single  vector  equation. 


dq 

Tt 


+  W-VJ  q  = 


-Vp 

P 


(13) 


We  now  substitute  for  q  with  the  gradient  of  the  velocity  potential  <D  as  defined 
in  equation  (8).  In  addition,  equations  (10),  (11)  and  (12)  are  used  to  obtain 
equation  (14)  from  equation  (13). 
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|[VO]+v 


-Vp 

p 


(14) 


It  is  desirable  to  express  the  right  hand  side  of  equation  (14)  purely  in  terms  of  a 
gradient  operator.  As  an  intermediate  step  we  may  see  that 


Yi 

p 


(15) 


where  A  is  the  dummy  or  umbral  variable  of  integration.  If  this  is  not  clear,  the 
explanation  follows  here.  From  equation  (5)  we  see  that  density  p  can  be 
evaluated  as  a  function  of  pressure  alone  (independendy  of  x,  y,  z  or  t)  such  that 
p  =  p(p)  explicidy.  The  lower  limit  of  integration,  pQ  ( t )  is  not  spatially 
variable.  The  upper  limit  is  p  (x}  y,  z,  t ) .  How  do  we  prove  that  equation  (15)  is 
correct?  This  is  shown  by  using  Leibnitz’s  Rule1  on  the  right  hand  side  of 
equation  (15),  to  obtain  the  left  hand  side.  Equation  (15)  is  often  written  in  a  less 
rigorous  form 


T  =  VJf  (16) 

Substitute  equation  (16)  into  the  right  hand  side  of  equation  (14).  Next,  reverse 
the  order  of  integration  with  respect  to  time  and  space  on  the  left  hand  side  of 
equation  (14).  Combining  all  the  terms  gives 


V 


d<X> 

It 


=  o 


(17) 


1.  Ltibcitt’*  Rule  ii  given  here.  Keep  in  mind  that  X.  if  m  umbrel  v arable  and  will  not  be  treated  u  a  function  of  x 
See  page  (365)  of  Wideband. 

s  /«**>*-  f 

A(*)  A  (a) 
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We  can  interpret  equation  (17)  as  a  single  vector  equation  with  three  components 
or  as  three  scalar  equations.  The  differentiated  quantity  in  each  of  the  three  sca¬ 
lar  equations  is  identical.  Clearly,  the  three  derivatives  are  zero  and  therefore,  the 
differentiated  quantity  is  independent  of  x,  y  and  z.  We  can  obtain  a  single 
expression  which  states  this  more  direcdy.  Integrate  each  of  the  three  vector 
components  of  equation  (17)  with  respect  to  x,  y  and  z  independently.  In  each 
case,  a  constant  of  integration  is  added  which  is  independent  of  x,  y  and  z 
respectively.  Since  the  same  quantity  must  result  from  the  integration  of  each  of 
the  three  vector  components,  the  constant  of  integration  must  be  independent  of 
x,  y  and  z.  The  only  variable  left  is  t  and  the  constant  of  integration  is  a  function 
of  time  alone.  Thus,  the  integration  of  (17)  leads  to  the  well  known  Kelvin’s  (the 
unsteady  version  of  Bernoulli’s)  equation. 

?  +  <18> 

We  are  left  with  the  task  of  determining  the  meaning  of  the  function  F  (/)  which 
arose  as  a  result  of  mathematical  manipulation  and  without  physical  insight 
Now,  equation  (18)  must  be  applicable  to  the  whole  flow  field  and  to  any  part  of 
it  We  now  specify  a  far  field  condition  (the  region  far  from  the  disturbance). 
Here,  the  flow  is  steady  and  the  streamlines  are  straight  Thus,  <l>  is  time 
invariant,  the  pressure  is  constant  and  the  velocity  is  assigned  a  constant 
magnitude  of  U .  So  our  far  field  conditions  are  q  =  U  ,  dp  =  0  and  the 
derivative  =  0.  By  restricting  equation  (18)  to  these  far  field  conditions,  we 
obtain 


F  (f)  =  C/2/ 2  (19) 

We  can  redefine  the  velocity  potential  such  that 

t  2 

♦  =  =  <b-~  (20) 

0 
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This  has  no  effect  in  the  interpretation  of  the  velocity  vector  and  the  substitution 
is  made  in  equation  (18).  The  resulting  equation  is  the  modified  Bernoulli’s 
equation. 


(21) 


We  momentarily  put  equation  (21)  aside  and  consider  the  continuity  equation. 
The  continuity  equation  (1)  can  be  put  into  vector  form 


dp 

37 


+  (<7-  Vp)  +p(V-  q)  =  0 


(22) 


Divide  by  p  and  substitute  q  -  V<J>  from  equation  (8)  and  equation  (20). 


i  r  dp 
p[37 


(V<t»Vp) 


+  (V*  v<i>)  =  o 


(23) 


We  now  step  back  and  see  what  we  have  accomplished.  Euler’s  formulation  has 
been  reduced  from  five  equations  and  five  unknowns  to  a  system  of  three  equa¬ 
tions  and  three  unknowns.  The  three  equations  are,  (23),  (21)  and  (5).  The  three 
unknown  system  variables  are  <fi,  p  and  p.  While  the  three  velocity  components 
no  longer  appear,  they  can  be  obtained  from  the  potential  4>.  We  still  have  the 
goal  of  formulating  a  single  equation  in  terms  of  a  single  system  variable  and 
independent  of  the  others.  We  will  choose  the  velocity  potential  0  as  this  single 
system  variable.  The  other  variables,  p  and  p,  can  be  obtained  subsequendy  to 
tire  solution  for  $• 

In  attaining  this  single  equation,  we  start  with  equation  (23)  which  has  three 
parts.  The  third  part  is  already  a  function  of  0.  The  first  part  can  be  made  a  func¬ 
tion  of  <f>  with  the  following  manipulations  based  on  equations  (5)  and  (6)  and 
Leibnitz’s  Rule. 
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Interchanging  the  first  and  last  steps  gives  a  useful  intermediate  formula. 


a  f  dp  ^ 
3rJ  p 


(24) 


The  left  hand  side  of  equation  (24)  employs  abbreviated  notation  which  is 
consistent  with  the  meaning  used  earlier  in  equation  (16).  We  see  where  the  right 
hand  side  of  equation  (24)  fits  in  a  time  differentiated  form  of  equation  (21). 
Making  the  substitution  in  equation  (21)  we  obtain 


a 

~9 


5p 

Tt 


d 

Tt 


q 

Tt*l 


2l 


Divide  by  a  to  obtain 


irap- 


r-i 


a 


a 


a<j)  o 

Tt+2 


2-1 


(25) 


Equation  (25)  can  be  substituted  into  the  left  hand  side  of  equation  (23).  Finally, 
the  second  term  of  equation  (23)  can  also  be  put  in  terms  of  <J>  with  the  following 
manipulations.  Starting  with  equation  (21)  we  obtain 


-V 


a<j> 

Tt 


(26) 


First,  using  equation  (16)  and  then  equations  (5)  and  (6),  we  see  that  equation 
(26)  can  also  be  written  as 


-V 


^  ,  <t 
Tt  +  i 


2 1 


=  Ye  =  l^Evp  =  ^vP 

p  p<*p  p 


Interchanging  the  first  and  last  step  and  dividing  by  a2  gives 


Ye  - 

p 


-i 

a2 


3$  7 
Ti  +  2 


2-1 


(27) 
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Finally,  taking  the  dot  product  with  V<|»,  it  follows  that 


(V<j>.  Vp) 


f — 1 
2 

a  j 


V(J>  *  V 


** . « 

dt  2 


2 1 


(28) 


This  is  in  the  desired  form,  ready  to  substitute  into  equation  (23).  But,  first  we 
cany  out  the  operations  on  the  right  hand  side  of  equation  (28)  to  give 


(V^Vp)  -± 
P  *  a 


d  2 

V<|>  ^(V<j>)+V4>  V(|) 


Now,  use  equations  (8)  and  (9)  and  note  that  ( q 1  =  V<J>  •  V<|>) 


(V<{)  •  Vp)  _i 


a2  l 


d\q2 

dt 


+  V<|>- V 


211 


q_ 

2 


(29) 


Equations  (25)  and  (29)  can  be  substituted  into  equation  (23)  to  give  the  full 
potential  equation. 


vV-i 

a 2 


§+|^+v*'v 


V 

2 


=  0 


(30) 


Equation  (30)  is  almost  in  the  desired  foim  of  a  single  equation  in  $  alone.  We 
address  the  parameter,  a  -  a(x, y,  z,  t )  in  die  following  development  Starting 
with  equation  (21),  we  substitute  for  p  (p) . 


d<J> 

Tt 


Y-l 


tp*-1’] 


(31) 


Using  equations  (5)  and  (6),  we  see  a2  =  yp^~  ^  and  substitution  in  equation 
(31)  gives  the  following  result 
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% 


9<f> 

Tt 


V  Y-l  J 


(32) 


Here,  a0(t)  is  the  speed  of  sound  associated  with  the  reference  (e.g.:  the  far 
field)  conditions  p  =  p0  and  p  =  po.  So,  by  equation  (32)  we  have  a  formula 
for  a  in  terms  of  <f>.  This  formula  for  a  can  be  used  in  equation  (30)  to  obtain  the 
one  single  equation  in  <}>.  However,  the  form  of  this  equation  is  complicated  and 
we  don’t  solve  it  anyway.  At  this  point,  we  assume  a  is  somehow  restricted  and 
address  this  issue  later.  We  expand  equation  (30)  using  subscripts  to  denote 
differentiation. 


v20  - 


la 


(t>  - 


r! 

a1  j 


r  l -j 
la  J 

[- n 

■y 


+  +  fk) 


(33) 


The  steady  state  form  of  equation  (33)  is  obtained  by  setting  to  zero  the  deriva¬ 
tives  with  respect  to  /.  We  proceed  under  the  assumption  that  a  steady  state  solu¬ 
tion  to  this  non-linear  equation  exists. 


We  now  turn  our  attention  back  to  equation  (33).  The  unknown  <j>.  is  divided  into 
two  components,  a  steady  state  component  (bar)  and  a  small  disturbance 
component  (tilde)  which  is  time  dependent  Likewise,  p  and  p  will  be  divided 
into  a  steady  state  component  and  a  small  disturbance  component  (tilde). 


$  (x,  y,  z,t)  =  $  (x,  yf  z)  +  $  (*,  y,  z,  t ) 

(34) 

p  (x,  y,  z,t)  =  p(xty,z)  +p  (x,  y,  z,  t ) 

(35) 

p  (r,y,  2,  /)  =  p(*,y,r)  +  p(x,y,z,f) 

(36) 

The  speed  of  sound  is  assumed  time  invariant  in  this  linearization  process.  So  we 
denote  this  restriction  as 
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a  (x,  y,  z,  t)  =  a  (x,y,  r)  (37) 

(Ultimately,  we  will  keep  a  constant.  One  will  see  that  a  higher  order  approxi¬ 
mation  has  no  effect  on  the  resulting  linear  formulation.)  We  substitute  equations 
(34)  and  (37)  into  equation  (33)  and  delete  any  non-linear  terms  in  <p  and  its 
derivatives.  Furthermore,  we  subtract  out  the  steady  state  condition.  We  obtain 
the  following  time  linear  partial  differential  equation  (PDE). 

($XX  +  $yy  +  U  ”  I" ~2~|  +  2^A/ +  2  Vvr  +  A/  + 

La  J 

♦X  +  $X  +  $X  +  2  <WX  +  +  + 


2  +  M>xA  +  + 

2($,$X  +  *X^/  +  ^X)  ^  =  0  (38) 

Here,  we  can  clarify  our  restriction  on  the  speed  of  sound  a.  If  we  had  assumed 
that  a  was  of  a  higher  order,  the  high  order  terms  would  have,  dropped  out  when 
we  dropped  all  the  non-linear  terms. 

While  partial  differential  equation  (38)  is  linear  in  the  solution  for  is  still 
difficult  to  analyze  for  any  general  description  of  the  steady  state  field 
v  (x,  y,  z)  .  This  spatially  variable  4>  (*,  >»  z)  results  in  spaiiaiiy  variable  coeffi¬ 
cients  in  this  PDE  for  Certainly,  there  are  no  elementary  solutions  available 
for  the  entire  flow  field  described  in  equation  (38).  So,  we  choose  to  further 
restrict  our  PDE  to  simple  steady  mean  flows.  If  we  let  the  steady  mean  flow 
have  a  uniform  velocity  of  U  with  streamlines  in  the  x  direction,  then  for  the 
entire  flow  field,  the  coefficients  are  simply 

4>(.*,y,2)  =  Ux  (39) 

a(xty,z)  =  a0  (40) 

where  aQ  is  the  constant  value  of  a  in  the  far  field. 
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Substituting  steady-state  components  given  by  equations  (39)  and  (40)  into 
equation  (38),  we  obtain 


(4> 


XX 


aoJ 


n 


(♦«+*  utx'+r+j 


(41) 


By  collecting  similar  terms  we  produce  the  classical  linear  small  disturbance 
velocity  potential  PDE.  We  now  name  this  equation  the  aerodynamic  potential 
equation  given  here  as  equation  (42).  We  have  defined  the  Mach  number  as 
M  -  U/a.  Both  the  steady  velocity  U  and  steady  speed  of  sound  a  were 
assumed  to  be  constant  throughout  the  flow  field  in  the  process  of  arriving  at 
linear  equation  (41).  For  a  non-linear  solution,  these  quantities  would  not  be 
constant  It  follows  that  for  linear  small  disturbance  theory,  the  Mach  number  is 
assumed  constant  for  the  entire  flow  field.  Later,  we  will  use  die  notation, 
p2  =  1-Af2. 


* 

f 1  ^ 

W 

2 

Kao  y 

(42) 


Now,  all  the  coefficients  are  constant  and  we  can  identify  elementary  solutions. 
This  will  be  assumed  to  be  the  governing  PDE  for  describing  the  aerodynamic 
behavior.  The  boundary  value  problem  is  comprised  of  equation  (42)  ana  the 
linear  boundary  conditions  to  be  specified  in  Section  IV.  It  can  be  shown  that  the 
solution  to  this  boundary  value  problem  is  unique. 
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The  Linearized  Pressure  Equation 


We  can  solve  for  <}>  (x,  y,  z,  t )  using  the  boundary  value  problem  comprised  of 
the  aerodynamic  potential  equation  (42)  and  the  boundary  conditions  to  be 
described  later.  However,  we  are  really  interested  in  the  pressure.  For  this 
reason,  we  use  Bernoulli’s  equation  (21)  to  develop  a  linear  expression  for  p  as  a 
function  of  <j>.  In  other  words,  having  developed  a  boundary  value  problem  for 
$  (x,  y,  r,  /) ,  independent  of  the  pressure  p,  we  now  develop  a  linear  formula  for 
determining  pressure.  This  functional  relationship  will  be  given  as  equation  (52). 
Furthermore,  we  will  express  $  as  a  function  of  p  in  equation  (73). 

Using  equation  (5)  we  can  write  pressure  as  an  explicit  function  of  density. 


Thus  the  integral  expression  in  equation  (21)  can  be  written 

p 

1 

p. 

Carrying  out  the  integration  on  the  right  hand  side  of  equation  (44),  we  obtain 


dp  _  yp0 


VPotv-2. 

-/f-JP  dfi 


(44) 


The  Linearized  Pressure  Equation 


Equation  (45)  can  be  substituted  into  integral  equpuon  (21)  to  obtain  the  non¬ 
linear  algebraic  expression  for  pres.- ore  in  terms  of  <j>  alone.  In  equations  (35) 
and  (36)  we  introduced  the  small  disturbance  notation  for  pressure  and  density. 
Here,  we  will  linearize  about  the  constant  far  field  pressure  condition 
p  (x,  y,  z)  =  p0  and  density  p  (x,  y,  z )  =  p0.  Since  we  are  interested  in  linear 
aerodynamics,  we  use  the  linear  part  of  a  Taylor  equation  identified  here  for 
some  function  F  (p)  synonymous  with  the  right  hand  side  of  equation  (45). 


p 

j~  =  F  (P)  *  F  iP0)  +  F'  (p0)  ip -p0)  +  hot  (46) 

Carrying  out  these  operafi  is  on  equation  (45)  and  simplifying  gives 


iP-Po) 


(47) 


This  linear  expression  can  be  substituted  into  equation  (21).  The  derivative  <j>;  in 
equation  (21)  is  easily  linearized  as 


90  90 

Ttadi 


(48) 


2 

Finally,  the  term  q  / 2  in  equation  (21)  is  linearized  about  the  free  stream 
velocity  components  u  =  U,  v  =  0,  and  w  =  0.  Starting  with 


2  i 

q  1,2,2,  2v 
2  =  21'"  +W) 


(49) 


We  use  a  multi-variable  Taylor  series  expansion  on  equation  (49)  to  obtain 


^*X-U2  +  U(u-U)  =  + 


(50) 
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We  substitute  equation  (47),  (48)  and  (50)  into  Bernoulli’s  equation  (21)  to 
obtain  the  linearized  expression 


(P-Po) 


(51) 


It  may  appear  that  this  expression  is  not  satisfied  at  the  far  field  condition  where 
<j>f  w  0,  =  0  and  p  =  p0.  This  is  easily  explained.  In  the  development  of 

equation  (21),  we  defined  <j>  as  an  alternative  definition  of  the  velocity  potential 
O.  In  doing  so,  we  essentially  set  the  far  field  pressure  to  zero.  If  p0  is  to 
describe  the  actual  far  field  pressure,  then  we  subtract  the  constant  p  (U2/ 2) 
factor  and  we  obtain  the  result 


(P-Po)  =  ~Po  <♦, +  UV  (52) 

It  turns  out  the  p^C/  /2  factor  will  be  inconsequential  because  we  will  be 
primarily  interested  in  die  pressure  differences  between  the  upper  and  lower 
surface  of  a  wing. 

Equation  (52)  is  a  very  important  formula.  In  Section  n,  a  linear  boundary  value 
formulation  for  the  potential  (j>  (x,  yt  r,  r)  was  derived  from  Euler’s  equations. 
Once  the  potential  function  <j>  (x,  y,  z,  t)  is  determined,  equation  (52)  is  used  to 
determine  the  pressure  p  (x,  y,  z,  t) .  Later  in  this  report,  it  will  be  useful  to  have 
a  formula  for  the  potential  given  the  pressure.  The  following  development1 
achieves  this  formula  as  equation  (73). 

We  start  with  equation  (52).  We  temporarily2  define  p  -  (p0 -p)  /po. 

p  =  $,  +  U$x  (53) 

We  use  the  method  of  characteristics  to  solve  this  problem.  This  method  utilizes 
a  coordinate  transformation  such  that  equation  (53)  becomes  an  ordinary  differ¬ 
ential  equation  which  can  be  easily  integrated.  The  new  coordinates  are 


1.  See  Ute  rcfemwe  by  Colioo. 

2.  P  deerfy  auepceuuie.Thif  previaw  deflation  i*  mi  to  be  coofcwd  with  Oi*  new  tea- 

paniy  um  of  Uk  tymbcJ  p. 
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£>-£>(*,  t) 

11  =  T1  (X,  t) 


(54) 

(55) 


Tlie  potential  function  and  the  pressure  function  are  given  new  designations  in 
the  new  coordinate  frame.  (Do  not  confuse  \y  here  with  die  acceleration 
potential  of  Section  XI.  No  connection  is  intended.) 

4>(jr(^,r|),r(^,  rO)  =\|>(4,ti)  (56) 

P  (X,  t)  -  p  (£,  n)  (57) 

We  ignore  the  role  of  the  y  and  z  coordinates.  It  will  be  seen  that  they  have  no 
influence  on  the  final  solution.  By  the  chain  rule,  we  have  from  equation  (56). 

♦,  =  Y&+Vnn,  (58) 

=  V&  +  VA  (59) 

Substituting  equations  (57),  (58)  and  (59)  into  equation  (53),  we  obtain 

P  =  (V&  +  MW  +y  (V&  +  V,n,)  (60) 

Rearranging  the  terms  gives 


f>  =  +  +  Vt|('l,+  (61) 

We  are  free  to  choose  the  relationship  between  the  (^,  Tj)  coordinates  and  the 
( x ,  y)  coordinates.  So  we  now  specify 


(il, +  (/’!,)  =  0 

(62) 

This  is  satisfied  if  we  simply  choose 

T)  (*,  t )  —  x—  Ut 

(63) 

Equation  (61  j  becomes 

i>  = 

(64) 

If  we  specify  ^  such  that 
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4  =  x  (65) 

we  obtain  the  relationship  between  pressure  and  the  potential. 

P  =  (66) 

Thus,  equation  (61)  has  been  reduced  to  equation  (66).  Now  we  simply  integrate 
equation  (66).  The  constant  of  integration  is  chosen  at  £  =  x  =  -<*>,  the  far 
upstream  condition. 

ft 

y(S,T|)  (67) 

— oo 


Now  we  work  to  cast  this  expression  in  the  original  coordinates.  From  equations 
(63)  and  (65)  we  have  the  inverse  relations 


x  =  5 


We  make  the  change  from  p  to  p  in  equation  (67). 


(68) 

(69) 


Now 

♦  Cm)  =  ¥(5*11)  =  W(x,x-Ut)  (71) 

Therefore,  by  substituting  equation  (71)  into  equation  (70),  we  obtain 

»«.*)  =  Jj  ]dX  (72) 

By  replacing  p  -  (pQ-p)  /p0  and  rearranging  terms,  we  obtain  the  desired 
result 
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(73) 


In  this  formula,  the  reader  is  reminded  that  X  is  die  dummy  variable  of 
integration  representing  integration  in  the  x  direction.  We  dropped  the  y  and  z 
dependence  as  a  matter  of  convenience.  Adding  y  and  z  back  to  the  argument 
list  of  the  pressure  function  p  (x,  t)  =*p  (x,  y,  z ,  t)  in  the  integrand  will  not 
change  the  integral  formula.  In  other  words,  equation  (73)  for  (*,y.  z,  t)  is 
written  as 


4>(*,y,  2*  0 


dX 


Again,  this  formula  is  used  if  we  know  p  (x,  y,  z,  t)  and  we  need  the 
corresponding  function  <j>  ( x ,  y,  z,  t) .  Now  X  appears  twice  in  the  argument  list, 
once  in  the  x  and  once  in  the  t  place  holders  of  the  pressure  function 
p  ( x ,  y,  z,  t) .  This  has  the  influence  that  as  we  integrate  downstream  from 
x  =  -oo,  we  evaluate  the  pressure  at  time  prior  to  t .  Later  in  this  report,  we  will 
refer  to  this  earlier  time  in  the  integrand  as  retarded  time.  Therefore,  if  we  know 
the  pressure  for  all  time  prior  to  time  r,  we  can  evaluate  the  above  integral 
expression  for  the  potential  <J>  (*,  y,  z,  t) . 
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Linearized  Boundary  Conditions  from  First  Principles 


The  aerodynamic  potential  equation  (42)  was  formulated  in  section  (II).  It  can  be 
shown  that  the  solution  to  this  partial  differential  equation  is  unique  given  the 
appropriate  boundary  conditions.  The  boundary  condition  specifies  the  potential 
or  a  directional  derivative  of  the  potential  on  all  surfaces  which  define  the 
computational  domain.  The  directional  derivative  of  die  velocity  potential  is  a 
component  of  the  velocity  vector.  So  we  specify  either  the  potential  or  a 
component  of  the  velocity  on  the  surface  of  the  computational  domain. 

When  the  flow  over  a  flight  vehicle  is  addressed,  the  computational  domain  is 
defined  interiorly  by  the  surface  of  the  flight  vehicle  and  the  trailing  wake  and 
exteriorly  by  the  far  field  conditions,  line  domain  may  be  nominally  fixed  with 
respect  to  the  vehicle  body  of  interest.  Certainly,  this  is  usually  the  case  for  most 
aerodynamic  developments  and  this  has  been  the  case  in  this  text  However,  one 
may  have  a  reason  to  attach  the  frame  of  reference  to  the  atmosphere  and  let  the 
vehicle  pass  through  the  reference  frame.  For  this  special  case,  the  interior  sur¬ 
face  of  the  computational  domain  moves  and  therefore,  the  boundary  condition 
is  applied  over  a  moving  surface.  Here,  we  will  only  address  boundary  condi¬ 
tions  on  a  surface  which  is  nominally  fixed  with  respect  to  the  frame  of  refer¬ 
ence. 

Vehicle  Surface  Boundary  Condition 

A  time  variant  surface  in  three  space  can  be  described  by  the  equation 

F(xty,z,t)  =  0  (74) 
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Of  course,  it  is  unrealistic  to  formulate  a  closed  fonn  expression  for  the  function 
F  if  the  surface  is  detailed  as  may  be  the  case  for  a  flight  vehic’e.  However,  one 
can  always  formulate  the  function  F  for  a  sufficiently  small  patch  of  the  total 
surface. 

The  boundary  condition  on  the  vehicle  specifies  that  the  flow  is  tangential  to  the 
surface.  In  other  words,  there  is  no  component  of  flow  normal  to  the  surface. 
Mathematically,  this  is  described  by1 

+  VF  =  0  (75) 

This  equation  can  be  linearized  about  any  reference  shape.  We  have  linearized 
the  potential  equation  (33)  about  an  undisturbed  uniform  flow  as  described  by 
equation  (39).  Thus,  the  boundary  condition  will  be  linearized  in  kind,  about  an 
undisturbed  and  uniform  flow.  As  mentioned  earlier,  this  is  a  severe  restriction. 
Basically,  this  limits  us  to  modelling  flow  disturbances  over  slender  bodies  and 
thin  wings. 

Here,  we  will  linearize  equation  (75)  specifically  for  a  thin  wing.  We  denote  the 
functional  description  of  the  surface  of  the  wing  as  F  =  Fw  (x,  y,  z,  t) .  This 
function  is  now  constrained  to  two  uncoupled  components,  the  deformation  of 
the  midplane  hm  and  the  thickness  envelope  ht  about  the  undeformed  midplane. 
The  undeformed  midplane  is  conveniently  designated  as  the  z  -  0  plane.  This 
is  stated  mathematically  as 

Fw  (x,  y,  z,  t)  *  z  -  hm  (x,  y,  t)  ±  ht  (x,  y,  t)  =0  ^ 

For  the  linearized  flow  about  a  uniform  free  stream,  V  =  Ui,  the  aerodynamic 
velocity  vector  is  mathematically  described  as 

V  =  ( U+u)i+  (v)j+  (w)k  (77) 


1.  Kirwacheti,  pp  191 
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For  the  remainder  of  this  text,  we  redefine  u,  v,  and  w  differently  than  in  equa¬ 
tion  (7)  to  represent  the  small  disturbance  from  the  uniform  free  stream.  We  sub¬ 
stitute  equations  (76)  and  (77)  into  equation  (75)  and  denote  h  -  hm±  hr 

-  ^  -  (U  +  u)  ^  +w  s=  0  (78) 

We  desire  a  linear  relation  between  the  velocity  components  at  the  surface  of  the 
wing  and  the  function  h  (.x,  y,  t) .  So,  we  make  equation  (78)  linear  in  h,  u,  v 
and  w  (keeping  in  mind  that  the  derivative  is  a  linear  operator)  by  dropping  the 
non-linear  terms.  Equation  (78)  now  becomes 


Bh 


+  U 


Bh 

Tx 


(79) 


It  is  now  clear  that  the  components  of  (h  =  hm  ±  ht)  can  be  treated  indepen¬ 
dently  in  the  linearized  boundary  conditions.  This  is  especially  important  when 
the  dynamic  response  of  a  wing  is  considered.  Here,  we  normally  assume  the 
thickness  effects  are  not  time  dependent. 


h  (x,  y,  t )  =  hm  (x,  y,  t )  ±  ht  (x,  y)  (80) 

Therefore,  when  analyzing  the  dynamic  response  of  a  wing,  we  superimpose  the 
dynamic  response  due  to  hm  on  a  separate  time  invariant  solution  using  ht.  This 
has  an  important  influence  in  our  choice  of  the  doublet  sheet  to  model  the  aero¬ 
dynamics  for  the  time  dependent  flow  over  the  wing  midplane  and  our  choice  of 
source  panels  to  model  the  aerodynamics  for  the  time  invariant  component  of 
flow  the  wing  thickness  envelope.  The  concepts  of  a  doublet  sheet  and  a  doublet 
lattice  will  be  discussed  later. 


The  Far  Field  Boundary  Condition 

The  far  field  boundary  condition  is  enforced  at  far  distances  from  the  interior 
boundary  where  the  flow  is  uniform.  It  will  be  clear  that  the  far  field  condition  is 
satisfied  automatically  when  one  uses  a  superposition  of  sources  or  doublets  on 
the  interior  boundary.  The  influence  of  sources  and  doublets  dies  out  at  infinite 
distances. 
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The  Trailing  Wake  Boundary  Condition 

Steady  state  lift  cannot  be  sustained  where  there  is  no  viscosity.  However,  the 
aerodynamic  potential  equation  (42)  is  restricted  to  irrotational  and  therefore 
inviscid  flow.  This  contradiction  is  corrected  by  fixing  the  flow  circulation  about 
the  airfoil  to  meet  the  Kutta  condition.  The  Kutta  condition  specifies  a  smooth 
and  finite  flow  off  the  sharp  tailing  edge  of  a  lifting  surface  in  incompressible 
flow.  The  velocity  vector  is  not  allowed  to  deflect  as  the  flow  passes  over  the 
trailing  edge.  If  it  does  deflect,  the  velocity  becomes  locally  infinite.  This  trailing 
edge  condition  and  wake  are  completely  characterized  for  incompressible  flow1. 

Linearized  steady  compressible  flow  over  planar  wings  can  be  transformed  to 
the  incompresible  case  (using  the  Prandtl-Glauert  transformation)  and  therefore 
the  trailing  edge  condition  is  well  understood.  For  linearized  unsteady  compress¬ 
ible  flow,  the  trailing  edge  condition  is  not  as  clearly  characterized.  For  instance, 
at  high  frequency,  it  is  experimentally  known  that  the  flow  off  the  trailing  edge  is 
not  tangential.  This  is  an  important  topic  and  warrents  further  study2.  For  linear¬ 
ized  flow,  we  introduce  a  trailing  wake  or  sheet  with  the  property  that  the  pres¬ 
sure  difference  across  the  sheet  is  zero.  This  is  the  only  condition  imposed  on  the 
wake.  We  allow  a  non-tangential  flow  on  the  wake.  A  few  of  the  consequences 
of  these  assumptions  are  now  discussed. 

For  our  linearized  boundary  value  problem,  the  wing  and  therefore  the  trailing 
edge  are  in  the  z  =  0  plane.  In  time,  the  mathematical  representation  of  the 
wing  slices  out  a  plane  region  as  it  moves  forward  through  the  air  at  velocity  U. 
We  will  treat  this  planar  region  as  the  trailing  wake.  (We  could  have  a  non-planar 
wake.  However,  for  our  linearized  system  of  equations,  the  added  accuracy  is 
not  warrented.)  It  is  assummed  that  there  is  no  jump  in  the  non-zero  velocity 
component,  w  ~  <}>z  across  the  wake.  Therefore,  <f>r  is  symmetric  with  respect  to 
z  in  its  first  approximation  (eg.  Fourier  Series)  with  respect  to  z  .  With  <j>z  sym¬ 
metric,  it  follows  that  <j>,  (j^  and  <j>f  are  all  antisymmetric  with  respect  to  z  in  the 
first  approximation.  Therefore  the  pressure 


1.  See  Sections  13-8  through  13-10  of  Karamcheti  for  a  discussion  of  the  wake.  Chapter  7  discusses  Sow  discontinuity 

2.  The  work  by  Guderley  partially  addresses  this  topic. 
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p  =  $t  +  U$x  (81) 

is  antisymmetric  with  respect  to  z  in  the  vicinity  of  the  wake.  Since  pressure 
jumps  are  not  admissible,  it  follows  from  antisymmetry  that  the  pressure  is  zero 
in  the  plane  of  the  wake. 

The  trailing  wake  boundary  condition  for  steady  flows  requires  <j)  to  be  antisym¬ 
metric  across  the  wake.  We  allow  the  possibility  of  an  anti-symmetric  jump  in  <j) 
across  the  wake.  While  continuity  in  pressure  is  a  requirement,  continuity  in  <j>  is 
not.  We  are  free  to  use  the  wake  as  a  boundary  on  the  domain  of  the  potential 
field.  Since  <}>  is  antisymmetric  and  discontinuous  across  the  wake  boundary,  so 

Boundary  Conditions  on  an  Oscillating  and  Deforming  Wing 

For  simple  wings,  it  is  often  expedient  to  represent  the  deformations  in  terms  of 
polynomials  in  space  and  harmonic  in  time.  Certainly,  this  is  the  case  for  flow 
over  a  plate.  We  assume  a  polynomial  of  order  nxinx  and  ny  in  y.  Frequency  is 
denote  as  ©  and  we  use  complex  notation.  So  we  constrain  the  out-of-plane 
deformation  to  the  following  series  expression: 


h  (x,  y,  t)  = 


nx  ny 


XXV/ 

j  =  o*  =  o 


teat 


(82) 


Of  course,  we  really  mean  to  equate  h  with  the  real  part  of  the  right  hand  side  of 
equation  (82).  We  substitute  equation  (82)  into  equation  (79)  to  obtain 


dh  rM 
w=T,+uTx 


=  10) 


nx  nv 


xxv/ 

7  =  o*  =  o 


+  U 


nx  »y 


XXW 


L/s=0*®0 


icat 

e  (83) 


J 
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The  complex  modulus  of  w  is  denoted  as  w  such  that  w  ~  weicat.  From 
equation  (83)  we  obtain 

w  =  X  '£aj*(i<s»Jyk  +  Vjx!~1yk) 

k  =  o  -j~  1 1  =  o 

The  main  reason  for  developing  equation  (84)  is  to  provide  an  example  of  the 
boundary  condition  formulation  which  may  be  used  in  the  doublet  lattice 
method.  Input  for  the  example  doublet  lattice  program  in  Appendix  A  is  in  this 
form. 
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Transformation  to  the  Acoustic  Potential  Equation 


In  Section  n,  we  started  with  the  Euler  equations  and  derived  the  aerodynamic 
potential  equation. 


(85) 


We  now  show  the  relationship  between  equation  (85)  and  the  classical  acoustic 
potential  equation. 


4 


4 


*4-  (h  +  <b 
-yj<,  -V, 


(86) 


It  turns  out  that  the  elementary  solutions  to  the  acoustic  potential  equation  are 
useful  in  identifying  other  elementary  solutions  to  the  aerodynamic  potential 
equation.  This  is  taken  up  in  the  next  section.  There  are  two  ways  in  which  one 
may  obtain  equation  (86)  from  equation  (85).  However,  the  interpretation 
differs. 


First,  we  note  that  in  the  derivation  of  equation  (85),  we  linearized  the  potential 
about  a  uniform  flow  with  velocity  U  in  the  positive  x  direction.  (Or,  what  is  the 
same,  the  body  moves  in  the  negative  x  direction.)  If  we  assume  the  flow  has 
zero  velocity,  then  equation  (85)  takes  the  form  of  equation  (86)  and  §  is 
identical  to  $  and  so  are  die  coordinate  frames. 
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However,  if  we  take  a  second  approach,  we  see  that  equation  (86)  can  be 
obtained  from  equation  (85)  by  the  follow.  \g  simple  translation  (also  known  as 
the  Gaussian  transformation). 


=  x-Ut 


y0  =  y 


zo  =  z 


T  =  t 


We  see  that  the  (x^0z0)  frame  moves  with  velocity  Ui  with  respect  to  the 
(x,y,z)  frame  and  is  therefore  motionless  with  respect  to  the  undisturbed 
atmosphere.  Next,  we  state  that  the  potential  in  the  (x,  y,  z)  frame  is  the  same  as 
the  potential  in  the  (x0,  y0,  z0)  frame.  We  distinguish  between  the  functional 
descriptions  with  an  over-tilde  and  an  under-tilde. 

( x ,  y,  z,t)  =0  ( x0 ,  y#  zQ ,  t)  (91) 

We  use  the  chain  rule  and  equations  (87)  through  (90)  to  cany  out  the 
differentiation  process.  We  denote  differentiation  with  subscripts. 

=  =  <92> 

4>,  =  ^  (*„)>  +  =  -(/$,.  +  4>t  (93) 

We  follow  through  with  less  detail  for  the  higher  derivatives  using  both  the 
product  rule  and  the  chain  rule  for  differentiation. 


A  =  A 
~xx  li 


♦«  =  -V*v.**v 

♦«  =  U'*V.-2U?1.X  +  ^ 

Of  course,  there  is  no  change  for  derivatives  in  the  lateral  directions. 
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(97) 


(98) 


By  substituting  equations  (94)  through  (98)  into  equation  (85),  we  obtain 
equation  (86). 


The  interpretation  of  the  two  approaches  is  different  The  first  approach  is  trivial. 
We  simply  set  the  velocity  to  zero.  The  second  approach  is  a  linear  transforma¬ 
tion.  The  aerodynamic  potential  equation  (85)  is  identically  the  same  as  the 
acoustic  potential  equation  in  a  uniformly  moving  frame. 
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The  Elementary  Solution  to  the  Acoustic  Potential  Equation 


in  the  previous  section,  we  showed  how  the  acoustic  potential  equation  (86) 
relates  to  the  aerodynamic  potential  equation.  In  this  section  we  seek  an  elemen¬ 
tary  solution  to  the  acoustic  potential  equation.  This  is  given  in  equation  (109). 
We  introduce  the  linear  Laplace  operator  V  and  restate  the  acoustic  potential 
equation  (86)  as  equation  (99).  (we  drop  the  tilde  underscore.) 


1  3^6  r-f2 

-2-i=  V2*  (99) 

a  dt 

In  equation  (99),  we  have  not  specified  <|>  in  any  particular  frame  of  reference. 
We  now  introduce  the  use  of  spherical  coordinates  where  r  is  the  radial  measure, 
0  is  the  measure  of  longitude  and  X  is  the  measure  of  latitude. 

<J>  =  <J>  (r,  8,  X,  t )  (100) 

The  acoustic  equation  (99)  takes  the  following  form1  in  spherical  coordinates. 


r  1  i 

3 

'sinJuf 

4. 

r  l  i 

_  r2sinX,_ 

5X 

T 

r2  ( sinA,) 1 

d2<|> 

ae2 


(101) 


We  are  looking  for  an  elementary  solution  to  the  acoustic  potential  equation 
This  elementary  solution  may  then  be  used  in  computing  complex  solutions  by 


1.  Hildebrand,  page  (313) 
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the  principle  of  superposition.  More  importantly,  the  elementary  solution  may  be 
adapted  to  solve  the  aerodynamic  potential  equation  (42).  We  make  an  educated 
guess  and  look  for  a  solution  to  equation  (101)  which  is  spherically  symmetric. 
We  temporarily  designate  this  solution  with  an  overbar  (not  to  be  confused  with 
the  steady  state  solution  in  Section  II). 


♦  M >(r,f)  (102) 

A  spherically  symmetric  flow  is  the  same  as  a  pulsating  source  (or  sink)  with 
radial  streamlines.  Physically,  such  a  flow  injects  mass  into  the  field.  The 
spherically  symmetric  form  of  the  acoustic  equation  is 


d2  <)> 

a? 


u2a 


Tr 


(103) 


In  order  to  maintain  the  undisturbed  far  field  condition,  we  seek  an  elementary 
solution  which  dies  off  as  r  — » <».  We  choose  the  lowest  order  expression. 


$(M)  (104) 

By  substituting  equation  (104)  into  equation  (103),  we  obtain  the  following 
hyperbolic  partial  differential  equation. 


(105) 


Hyperbolic  equations  have  characteristic  solutions1.  The  two  characteristic  solu¬ 
tions  to  equation  (105)  are  combined. 


/(r,  r)  =  /,(/■■*•  at)  +fe  (r  -  at)  (106) 


1.  *oe  page  399  of  Hildebrand 
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where  /■  and  fe  are  any  analytic  function.  Usually,  the  classical  acoustic  solution 
is  placed  in  the  following  form  which  is  obtained  from  equation  (106)  by  modi¬ 
fying  the  form  of  ft  and  fe. 


f(r,t)  =  s,|r+£ 


(107) 


One  can  easily  verify  that  equation  (107)  is  a  solution  of  equation  (105).  Now, 
substituting  equation  (107)  into  equation  (104),  we  obtain  the  result 


(108) 


If  we  plot  ge  as  a  function  of  time  (assuming  some  initial  pulse)  we  see  that  ge 
represents  an  expanding  wave.  Likewise,  gi  represents  an  imploding  wave.  We 
choose  gt  =  0  because  the  acoustic  phenomena  of  interest  here  takes  the  form 
of  an  expanding  wave  from  a  central  disturbance.  We  drop  the  subscript  e  and 
write 


<|>(r,r)  = 


(109) 


Y 

The  argument  (t- -)  is  called  retarded  time  because  it  accounts  for  the  delay 
between  the  time  the  radiating  disturbance  was  initiated  at  r  =  0  until  the  time 
it  reached  the  distance  r.  We  may  assume  that  the  function  g  ( t )  is  known. 


It  is  a  straightforward  matter  to  verify  that  equation  (109)  is  a  solution  to  the 
acoustic  potential  equation  (86)  given  in  the  previous  section  in  terms  of  carte¬ 
sian  coordinates.  First  we  establish  the  following  relations. 


,  2  2  2^ 
r  =  (x  +  y +z ) 


1/2 


X 

r  =  - 


(110) 

(111) 
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1  *■ 

_  1  X 

rxx  ~  ~Z  71 


(112) 


We  now  proceed  to  develop  expressions  for  the  derivatives  of  $  *  <j>  given  in 
equation  (109). 


1  = 


— 1 

l. 

.  r 

x  ‘ 

t 

r 

3 

Lr  J 

8 

t  — 
a  J 

2 

.ar  J 

8 

t  — 
a. 

(113) 


d>  = 

r  XX 


i 

>— t 

_ _ 1 

r 

3*2  1  1 

• 

r 

*2  ■ 

tt 

„  r 

5  ~  3 
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L  u  j 
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V 
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a  [ 
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8 
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L  a  J 
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8 
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(114) 


The  derivatives  with  respect  to  y  and  z  are  derived  similarly. 


<h  = 


3 1  i 

S  3 
r  r 


8 


r 

*~a} 


„  2  < 
jy  1 

4  2 

ar  ar 


8' 


r 

1  a] 


.2  i 


2  3 
a  r  J 


8' 


f- - 

a  j 


(115) 


<b  = 
Yzr 


3^_  J. 

r5  ~  r 


8 


t- - 

aj 


3  z2  1 


ar 


ar2  j 


8' 


t-l 
a  J 


2  -i 


2  3  , 
ar  j 


5' 


r-- 

aj 


(116) 


The  derivative  with  respect  to  t  is  directly  derived. 


K  -  7* 


(117) 


We  substitute  equations  (114),  (115),  (116)  and  (117)  into  the  equation  (86)  and 
we  see  that  it  is  satisfied.  Thus,  we  have  confirmed  that  =  $  is  an  elementary 
solution  to  the  acoustic  potential  equations. 

Now  we  come  to  an  important  point  Up  until  now,  we  have  assumed  the  point 
source  is  stationary.  If  the  source  is  moving  in  time,  relative  to  the  point  x0%  y0, 
z„,  then  the  definition  of  r  changes.  From  equation  (110),  we  now  have 
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o  a  1/2 

r  =  [(x(t)  -x0)  +  (y(t)  -y0)2+  ( z(t )  -z0)2]  (118) 

The  expression  for  <|>fr  in  equation  (117)  is  not  correct  for  the  moving  source  and 
the  acoustic  equation  (86)  is  no  longer  satisfied.  So  we  formulate  a  new 
elementary  solution  which  represents  a  moving  source.  This  is  accomplished  in 
the  next  section. 
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The  Moving  Source 


In  this  section,  we  obtain  the  elementary  formula  (145)  for  the  moving  source.  In 
the  next  section,  we  transform  this  solution  to  the  moving  frame  of  reference. 
Then  one  may  demonstrate  that  this  formula  does  indeed  satisfy  the  aerody¬ 
namic  potential  equation  (42). 

The  following  explanation  is  a  review  of  the  derivation  provided  by  Garrick1. 
We  superimpose  a  train  of  stationary  sources,  all  in  a  line.  These  sources  are 
pulsed  in  a  sequence,  thus  giving  the  same  effect  as  a  single  source  moving  at  a 
constant  speed.  One  may  think  of  this  as  a  motion  picture  film  of  a  moving 
source  derived  from  a  series  of  photographic  frames. 

The  pulse  function  6  (f)  is  defined  to  be  0  when  t*  0  and  to  be  1  when  t  -  0. 
We  can  define  the  pulse  function  more  eloquently  with  the  following  continuous 
function. 


f 

0 

6(f)  =  lim 

1 

1  -f  cos 

*  itf*r 

a  ->0 

2 

.  a  J. 

0 

t  <  -a 

-a£t£a 

t>a 


(119) 


The  pulse  function  has  the  following  effect  when  placed  in  the  integrand.  When 
we  view  the  integration  process  as  the  limit  of  a  series  summation,  the  pulse 


1.  S«ep*$c671  of  the  wodeby  Ganici 
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function  selects  only  one  value  of  the  integrand  and  reduces  the  other  series  ele¬ 
ments  to  zero. 


/(f)  =  J/(x)8(f-x)dx  (120) 

— oo 

Given  a  set  of  orthogonal  x0,  y0and  z0  axes,  fixed  to  the  atmosphere,  place  a 
series  of  stationary  sources  coincident  with  the  x0  axis.  Thir  is  shown  in  figure  1. 
These  six  sources  shown  are  located  at  x0  -  where  the  index  i  ranges  from 
1  <  /  £  6. 


(xo>y0’  zo) 


— a  c  > 

£3  ^  ^>1 

*0 


Figure  L  Acoustic  Sources  on  the  x0  Axis 

All  six  stationary  sources  are  assumed  to  act  independently.  Because  the  acoustic 
potential  equation  is  linear,  we  can  superimpose  the  elementary  formula  (109) 
relating  the  potential  at  (xQ,  yc,  z0)  due  to  a  source  located  at  x0  -  tod 
y0  =  =  0. 


^(x0,y0,  zo>  0 


(121) 
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where 


rf  =  (x0-%)2-Ky0)2+(zf  (122) 

If  g  is  the  pulse  function,  then  equation  (121)  represents  the  potential  which 
arises  due  to  the  six  sources  which  are  simultaneously  pulsed  at  t  -  0  with  unit 
strength. 


$(x0iy0,z0,  t) 


I 

7  =  1 


(123) 


In  equation  (123),  we  see  the  effect  of  the  retarded  time.  Even  though  all  six 
sources  are  pulsed  simultaneously,  the  effect  at  the  point  (x0,  y0,  z0)  is  picked 
up  (r/a)  later.  In  other  words,  six  simultaneous  pulses  are  transmitted  to  the 
point  (x0,y0,z0)  with  a  different  delay. 


Now,  instead  of  pulsing  simultaneously,  we  pulse  the  sources  in  a  sequence, 
starting  with  the  source  at  ^  and  ending  with  the  source  at  ,  Each  source  is 
pulsed  with  a  time  delay  of  t  -  relative  to  t  -  0.  Furthermore,  rather  than  a 
unit  strength,  each  source  is  pulsed  with  strength  Fi .  Equation  (123)  now  takes 
the  form 


♦  (  •Wo>-V')  =  X 


t  =  1 


LriJ 


(f,)8 


<‘-*-2. 


(124) 


Instead  of  six  point  sources,  we  may  have  many  point  sources  along  the  x0axis. 
In  fact,  we  may  define  a  continuum  of  sources  as  a  limiting  process  as  the  num¬ 
ber  of  point  sources  go  to  infinity.  (While,  this  argument  is  not  rigorously 
defended  here,  it  can  be  shown  that  the  desired  formula  relating  the  potential 
field  due  to  a  moving  source  does  satisfy  the  aerodynamic  potential  equation.) 
The  summation  process  of  equation  (124)  is  now  treated  as  an  integral 


47 


Ths  Moving  Source 


$(x0>y0>zo>t)  = 


r( x)  J 


F(t)8 


(t- t)  - 


r(T) 


<7 


<J?X 


(125) 


where  x  is  die  dummy  variable  of  integration  which  one  may  think  of  as  running 
in  parallel  with  time  t.  The  variable  x  points  to  some  place  in  time  for  which  a 
uniquely  tagged  point  (or  differential)  source  is  active.  We  devise  the  definition 
of  r  (x)  based  on  rf  given  in  equation  (122). 


r2(x)  =  (at0-^(t))2  +  ^+22  (126) 

F  (x)  dx  is  the  differential  strength  of  the  source  pulse  which  is  uniquely  identi¬ 
fied  by  the  time  delay  x.  Now,  let  the  sources  be  pulsed  at  a  uniform  rate  of  -U 
along  the  x  axis.  Furthermore,  we  specify  when  x  =  0,  there  is  a  pulse  at 
x  =  0.  In  other  words,  we  specify 


and  equation  (126)  becomes 


(127) 


r2(t)  =  (*0  +  11t)2  +  )>2  +  z2  (128) 

We  now  strive  to  evaluate  the  integral  equation  (125).  The  resulting  formula, 
equation  (145),  is  the  desired  elementary  solution  for  a  moving  source.  In  order 
to  take  direct  advantage  of  the  sifting  property  described  by  equation  (120),  we 
employ  the  following  change  of  variable  in  equation  (125)  in  order  to  simplify 
the  process. 


-0  =  f-x- 


r(x) 


Substituting  for  r  (x)  from  equation  (128),  we  obtain 


-Qs  ,--c-{±)aXo  +  Ur)2+yl  +  zl) 


1/2 


(129) 


(130) 
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dx  1 
d6 


d& 


(131) 


In  the  process  of  changing  the  variable  of  integration  from  x  to  0,  we  require 
expressions  for  x  (0)  and  dx/d&.  This  task  is  easier  than  one  may  suppose  at 
first  According  to  the  sifting  principle  of  equation  (120),  we  evaluate  the  inte- 
grand  at  ©  =  0  only.  This  gives  the  result 


4>  (*o»  yo>  zo>  0  = 


1  1 


L '•(*)_! 


F(x) 


dx 

le 


0  =  0 


(132) 


We  now  evaluate  equation  (132)  for  the  potential  which  arises  from  a  moving 
source.  This  is  achieved  in  equation  (145)  using  equations  (135)  and  (143). 

From  equation  (130),  obtain  the  following  quadratic  in  x  for  0  =  0. 


(|32)t2-2 

»2  i  i  -2 


Ux 

r+ 

2 

x  + 

L  a 

L 

K 


(*;;+>’»+ 4) 1 


a 


=  0 


(133) 


where  (3  =  \-M  and  M  =  U/a  is  the  Mach  number.  From  this  quadratic 
equation,  we  expect  two  solutions  for  x. 


t= 


^j±j((*.+c,()2+p^+f^)1/2 


(134) 


For  subsomc  flow  (Af  <  1) ,  we  choose  to  limit  t  <  t  which  limits  equation  (134) 
to  one  root1.  Equation  (134)  becomes 


1 .  Garrick  give*  to  excelled  graphical  grplaq.tj/i^  fa  jk, 
674  tad  676  respectively. 


mbaoaic  and  mipawxk  caaea  in  figu.net  F.4a  aod  F,4b  oc  ptgea 
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f  + 


where 


R  =  ((Xo  +  Ut)2  +  t>2yz0  +  tfz20)l/2 


(135) 


(136) 


Equation  (135)  is  the  expression  for  x  when  0  =  0  and  U<a.  We  now 
compute  dx/dB  with  ©  =  0.  Equation  (130)  is  rearranged,  squared  and 
differentiated  to  obtain  equation  ( 1 37). 


jjt[a2(T-,-e)2]  =  ^[(*0+t/t)2+^+zJ] 


Carry  out  the  differentiation  and  then  set  0  =  0. 


a  (x-r) 


dx 

dB 


~  1 


=  (xQ+Ux)U 


Now  solve  for  dx/dB, 


'  dx ' 

dB 


[a2(T-»)-tf(*„  +  E/t)]^  =  a2(t-») 

Next,  we  observe  from  equation  (129),  for  0  =  0 

r  =  a(t-x) 

Making  the  substitution  in  equation  (139)  gives 

[- ar-U(x0  +  Ux ))  =  -or 

1  dx  _  a 

rdB  ar  +  U(x0  +  Ux ) 


(137) 

(138) 

(139) 

(140) 

(141) 

(142) 
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Using  equations  (140)  and  (135),  we  substitute  for  r  and  x  on  the  right  hand  side 
of  equation  (142).  After  some  simple  algebraic  manipulation,  we  arrive  at  die 
following  simple  result 


1  __  1 
rdB  R 


(143) 


The  definition  for  R  was  given  in  equation  (136).  Finally,  we  use  equation  (143) 
in  equation  (132)  to  obtain 


$(x0’y0>zo>  0 
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R] 


9  =  0 


(144) 


From  equation  (135),  we  have  an  expression  for  x  when  0  =  0.  We  make  the 
substitution  in  equation  (144). 


$(xo>y<>>  zo>*) 
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where  R  is  defined  in  equation  (136). 

This  is  the  fundamental  solution  for  the  potential  which  arises  due  to  a  source 
moving  along  the  x  axis  with  constant  velocity  of  ~Ul.  In  the  next  section,  we 
transform  the  coordinates  to  a  moving  frame. 
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Tiie  Elementary  Solution  to  the  Aerodynamic  Potential  Equation 


Our  objective  of  the  past  three  sections  has  been  to  derive  elementary  solutions 
to  the  aerodynamic  potential  equation  (42)  which  may  be  used  to  model  the  flow 
over  wings  and  bodies.  In  Section  V,  we  recognized  that  the  aerodynamic  poten¬ 
tial  equation  is  related  to  the  acoustic  potential  equation  by  a  simple  Gaussian 
transformation.  Tire  coordinates  axes  of  the  acoustic  potential  equation  are  fixed 
to  the  atmosphere  while  the  coordinate  axes  of  the  aerodynamic  potential  equa¬ 
tion  move  with  constant  velocity  -Ul  relative  to  the  atmosphere.  The  elemen¬ 
tary  solution  to  the  acoustic  potential  equation  is  a  stationary  point  source  with  a 
spatial  decay  of  ( 1  /r) .  We  used  a  modified  form  of  this  solution  (109)  to  obtain 
an  elementary  solution  to  the  aerodynamic  potential  equation.  Then  a  complica¬ 
tion  arose.  We  discovered  that  a  simple  translation  of  the  stationary  source  in  the 
x  direction  does  not  satisfy  the  acoustic  potential  equation.  This  mathematical 
complication  is  the  result  of  compressibility  (also  referred  to  as  the  Doppler 
effect).  Here,  we  are  faced  with  the  apparent  compression  of  pressure  wave 
fronts  travelling  upstream  and  the  apparent  expansion  of  pressure  wave  fronts 
travelling  downstream.  So,  through  the  limiting  process  of  superimposing  a 
series  of  source  pulses,  we  simulated  a  constant  velocity  source  and  derived  the 
formula  (145)  for  the  resulting  potential. 

In  this  section,  we  apply  a  change  of  coordinates  to  the  moving  source  solution 
(145)  in  the  acoustic  frame  and  thereby  obtain  the  moving  source  solution  (152) 
in  the  original  constant  velocity  frame.  This  is  the  desired  elementary  solution  to 
the  aerodynamic  potential  equation.  One  may  directly  verify  that  equation  (152) 
solves  the  aerodynamic  potential  equation  by  direct  substitution. 


53 


Hie  Elementary  Solution  to  the  Aerodynamic  Potential  Equation 


Again,  our  objective  is  to  mathematically  model  the  flow  over  wings  and  bodies. 
One  approach  is  to  position  a  continuum  of  sources  on  the  wing  or  body  surface 
in  order  to  disturb  the  flow  and  thereby  satisfy  the  tangential  flow  boundary 
condition.  This  concept  of  a  spatial  continuum  of  sources  will  be  discussed  in 
die  next  section. 

The  moving  coordinate  frame  is  fixed  relative  to  the  structural  geometry  and  has 
a  velocity  -Ui  t dative  to  the  acoustic  coordinate  frame.  Therefore,  a  source 
moving  with  velocity  —Ui  in  the  acoustic  frame  is  now  fixed  relative  to  the 
moving  frame.  The  potential  which  arises  from  a  moving  source  was  presented 
as  equation  (145).  From  equation  (127)  we  know  that  at  t  =  0  the  single  moving 
source  is  located  at  the  origin  of  the  (xoy0zo)  axes.  We  may  modify  the 
elementary  solution  (145)  to  model  the  potential  4>  (x^y^)  due  to  a  single 
moving  source  which  maintains  a  constant  distance  (£,  q,  Q  relative  to  the 
source  located  along  the  x0  axis  at  (\o  =  -Ut) .  (We  revert  to  the  under-tilde  to 
denote  the  potential  in  the  stationary  acoustic  frame.  Furthermore,  the  functional 
notation  <j>  is  used  to  denote  the  transitional  state  between  the  stationary  and 
moving  frames.) 


^  y0>  *o'  4*  *1»  0  = 

and  from  equation  (136) 
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Now  we  use  equations  (87)  through  (91)  to  change  from  the  fixed  (xjoi0) 
coordinates  to  the  moving  (x,y,  z)  coordinates. 
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and 
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1/2 

R=  ((*-£)2  +  p2()>-tl)2  +  P2(*-4)2)  (149) 

We  make  algebraic  simplifications  to  equation  (148)  to  obtain  the  following 
equation. 


4*  (xo>  y o >  zo>  ,n>  o  — 


Ilf 


t  + 


M(jc-S) -R 


(150) 


As  a  final  step,  we  give  x  a  new  definition,  not  related  to  the  dummy  variable 
used  in  the  previous  section.  Here,  x  is  the  retarded  variable  and  it  represents  the 
time  delay  incurred  for  a  pulse  to  transit  from  its  origin  at  (§,  r\,  Q  to  the  point 
(x,  y,  z) . 


x  - 


-M  {x-%)  +  R  0 
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(151) 


The  form  of  equation  (150)  is  simplified. 


(*>y>  tri*  1 0 


U-x] 


(152) 


The  subscript  s  is  added  to  denote  the  source  solution.  Later  a  subscript  d  will 
denote  the  doublet  solution.  In  the  derivation  of  equations  (152),  (151)  and  (149) 
we  closely  followed  the  approach  taken  by  Garrick.  This  is  the  fundamental 
moving  point  source  solution  to  the  aerodynamic  potential  equation.  Thus,  our 
interim  objective  has  been  achieved.  Higher  order  solutions  can  (and  will)  be 
formulated  by  differentiating  equation  (152)  with  respect  to  x,  y  or  r. 

It  can  be  shown  by  direct  substitution  that  equation  (152)  solves  the  aerody¬ 
namic  potential  equation. 
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The  Source  Sheet 


In  the  previous  section,  we  showed  that  the  fundamental  source  solution  to  the 
aerodynamic  potential  equation 


2  U 

„2 


♦,r 
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=  0 


(154) 


is  given  by  the  following  simple  formula. 


<J>S  (xty,  z,  t) 


[t-x] 


(155) 


This  is  the  fonnula  for  the  potential  at  coordinates  (x>  y ,  z)  due  to  a  single  point 
source  at  coordinates  (£,,  r\,  Q  .  The  boundary  condition  for  the  flow  over  a  thin 
wing  was  given  in  equation  (79). 


dty  dh  dh 

M'  =  Tt  =Tt+u3i 


(156) 


where  h  (xt  y,  /)  describes  the  time  dependent  deformation  of  a  thin  wing  in 
the  (a,  y)  plane.  Tne  obvious  question  remains;  how  do  we  use  equation  (155) 
to  solve  for  the  flow  over  a  wing?  The  answer  is  not  simple  and  is  the  subject  of 
the  remainder  of  this  text.  We  still  need  to  formulate  the  source  doublet  in 
Section  X  and  then  we  formulate  the  pressure  doublet  in  Section  XI.  We  use  the 
concept  of  the  doublet  sheet  to  develop  ihe  integral  formula  in  Section  XU  In 
this  section,  we  are  introduced  to  the  concept  of  a  source  sheet  While  we  will 
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not  use  the  resulting  formulae  for  a  source  sheet,  the  concept  is  directly 
applicable  to  a  doublet  sheet 


If  we  restrict  our  problem  to  the  flow  over  a  steady  wing  with  no  deflection  (i.e. 
only  the  wing  thickness  is  a  consideration),  the  boundary  condition  equation 
(156)  simplifies  to 


_  tM 

dz  "  °dx 


(157) 


Equation  (155),  for  the  elementary  point  source  solution  with  temporally  con¬ 
stant  strength  A ,  simplifies  to 


<t>5  (x,  y,  z,  4,  T],  0 


r^nUOi 

R 


(158) 


Where  A  is  the  unknown  strength  of  the  source  at  the  coordinates  (!;,  tj,  Q .  In 
order  to  satisfy  this  boundary  condition,  we  may  superimpose  n  point  sources 
located  at  ( x ,  y,  z)  =  (£,,  rj.,  0) .  We  satisfy  the  boundary  condition  on  the 
wing  surface1  at  n  points  located  at  (Xp  yp  Zj) .  The  thickness  envelope  is  sym¬ 
metric  above  and  below  the  (x,  y)  plane  so  it  is  sufficient  to  satisfy  the  bound¬ 
ary  condition  on  the  top  surface  only.  We  differentiate  equation  (158)  with 
respect  to  z  and  substitute  the  result  into  equation  (157)  for  each  static  point 
source. 
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So  we  have  a  system  of  n  equations  and  n  unknowns  which  can  be  solved  with 
linear  algebra.  One  expects  the  accuracy  of  the  solution  to  increase  as  the 
number  of  point  sources  and  control  points  is  increased.  We  can  reformulate 
equation  (159)  as  an  integral  if  we  consider  the  source  in  a  limiting  process. 


1 .  By  Mlufying  the  boundary  condition  on  the  wing  lurface,  we  are  inconhaent  with  our  auumed  line arizalian  a  the  wing 
midplane.  However,  one  expect!  the  boundary  coodition  on  the  wing  turf  ace  to  be  more  accurate  than  the  midplaae.  Betide*, 
wc  avoid  problem!  with  npgularilie*. 
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S 

Here,  we  have  discretized  the  ( x ,  y)  plane  into  differential  areas  of  size  dS  and 
located  at  *  =  ^  anti  y  =  q.  Each  differential  area  has  a  source  strength  of 
AdS.  Again,  the  steady  function  h  is  evaluated  at  (x,  y) .  The  radial  measure  R 
is  now  defined  as 


R 


dS 


(160) 
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1/2 


(161) 


If  one  is  given  the  value  of  h  (x,  y)  at  m  points,  equation  (160)  can  be  approxi¬ 
mately  solved  for  the  unknown  function  A  (4,  T|)  if  A  (£,  r|)  is  defined  in  terms 
of  m  approximating  components  with  constant  coefficients.  For  instance,  we 
may  form  a  composite  function,  A  (4,  T|)  by  superimposing  m  polynomials  in  £ 
and  t|,  each  polynomial  weighted  by  a  constant  (but  not  yet  specified)  quantitiy. 
We  integrate  equation  (160)  for  each  polynomial.  This  results  in  a  linear  system 
of  m  equations  with  m  unspecified  constants.  Alternatively,  we  may  spatially 
discretize  the  wing  planform  and  approximate  A  (£,  r\)  with  a  continuous  spline 
function  with  m  unspecified  coefficients.  Again,  equation  (160)  is  integrated  to 
obtain  a  linear  system  of  m  equations  with  m  unknown  coefficients.  We  may 
think  of  this  aerodynamic  modei  as  a  linear  system  with  m  independent  inputs 
( h  (x,  y)  at  m  points  on  the  surface)  and  m  dependent  outputs  (m  polynomial 
constants).  Once  the  approximating  solution  for  A  (S,,  q)  is  obtained,  the  poten¬ 
tial  is  determined  using  equation  (158).  Then  we  use  the  time  invariant  term  of 
equation  (52)  to  solve  for  the  pressure. 


F(*,y)  -  -pU 
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062) 


Thus,  one  example  of  how  one  may  use  the  elementary  point  source  solution  to 
obtain  a  continuous  solution  has  been  given.  Keep  in  mind  that  we  are  solving 
the  small  disturbance  problem  and  that  the  solution  breaks  down  at  stagnation 
conditions. 
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It  turns  out,  that  a  single  source  sheet  cannot  generate  a  pressure  difference 
across  the  (x,  y )  plane.  Therefore,  no  lift  can  be  generated.  Mathematically,  this 
is  seen  when  one  recognizes  the  symmetry  of  the  potential  above  and  below  the 
(x,  y )  plane.  For  this  reason,  we  investigate  the  source  doublet  in  the  next  sec- 
tion. 
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The  Source  Doublet 


In  Section  IV,  the  tangential  flow  boundary  condition  over  a  thin  wing  was 
separated  into  two  linear  components,  wing  thickness  and  wing  deformation. 
Thus,  we  can  treat  the  linear  boundary  value  problem  for  a  thin  wing  as  two 
separate  problems,  the  potential  which  arises  due  to  the  thickness  envelope  and 
the  potential  which  arises  with  the  deformation  of  the  wing  midplane.  The  total 
solution  is  the  superposition  of  the  two  component  solutions. 

In  the  analysis  of  most  linear  systems,  one  considers  the  steady  state  condition 
and  then  superimposes  the  time  dependent  response.  The  steady  state  solution 
for  a  wing  is  a  superposition  of  the  pressure  due  to  thickness  and  the  pressure 
due  to  steady  deformation  of  the  midplane.  The  time  dependent  response  is  asso¬ 
ciated  with  the  time  dependent  deformations  of  the  wing  midplane  alone.  This  is 
an  important  consideration  because  we  can  mathematically  model  the  flow  over 
a  deforming  wing  with  an  infinitely  thin  sheet 

In  the  previous  section,  we  demonstrated  a  solution  technique  using  a  source 
sheet  However,  it  was  pointed  out  that  due  to  the  symmetric  nature  of  the  source 
sheet  (with  respect  to  the  (x,  y)  plane),  it  was  not  possible  to  develop  a  pressure 
differential.  Therefore  no  lift  can  be  generated  with  a  single  source  sheet 
However,  it  is  possible  to  develop  a  pressure  difference  if  two  source  sheets  are 
placed  in  parallel.  This  can  cause  numerical  problems  if  the  two  sheets  are 
brought  close  together.  This  is  not  to  say  tlus  has  not  been  done.  On  the  contrary, 
there  are  many  examples  where  this  is  exactly  what  is  done.  However,  for  our 
linear  analysis,  this  results  in  a  waste  of  computational  resources.  Instead,  we 
can  formulate  the  limiting  condition  as  two  source  sheets  with  opposing 
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strengths  are  brought  close  together.  This  is  the  source  doublet  sheet  3h  this 
Section,  we  describe  the  elementary  point  doublet  formula. 

Given  the  linear  aerodynamic  potential  equation 
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and  the  fundamental  source  solution  of  equation  (152). 

♦,  =  j =/('-*) 

we  show  that  the  elementary  solution 
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is  also  a  solution.  We  substitute  equation  (165)  into  equation  (163)  to  obtain 
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Next,  the  order  of  differentiation  is  changed 
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The  term  in  the  square  brackets  is  known  to  be  zero  from  equation  (163).  Equa¬ 
tion  (167)  reduces  to 


Tt 
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Thus,  we  have  shown  that  equation  (165)  is  also  a  solution  to  the  aerodynamic 
potential  equation.  We  call  this  the  source  doublet  or  just  the  doublet  solution.  In 
order  to  give  some  physical  significance,  we  investigate  the  source  doublet  for 
the  steady  incompressible  case  (set  M  =  0).  Here  we  have  the  point  source 
solution  for  a  source  of  unit  strength  located  at  the  origin.  The  potential  is 

♦,  =  1  (169) 

with 

2  2  .  2  .  2 

r  ~  x  +y  +z  (170) 

According  to  equation  (165),  the  potential  which  arises  from  a  point  doublet 
solution  can  be  obtained  by  differentiation  of  the  source  solution  with  respect  to 
z. 


(171) 


Now,  we  take  a  second  approach  to  arrive  at  equation  (171).  We  bring  two 
sources  together  from  above  and  below  the  z  =  0  plane  as  shown  in  figure  (2). 
The  strength  of  the  sources  are  opposite  and  inversely  proportional  to  the  dis¬ 
tance  between  them.  Using  equation  (169)  we  obtain  the  combined  potential 


lim 
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IX 


Equation  (170)  is  modified  for  each  source  as  follows. 


(172) 


r2  = /  +  /+ (z  +  Q2  (173) 

r\  =*2  +  A(z-02  (174) 

We  temporarily  assume  z  *  0.  We  can  see  that  as  £  ->  0,  we  have  a  zero  in  the 
numerator  and  a  zero  in  the  denominator. 
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We  now  take  the  limit  as  £  -»  0  to  obtain  equation  (179).  Furthermore,  we  now 
allow  z  to  go  to  zero.  We  see  that  equation  (179)  agrees  with  equation  (171). 

*d  =  ^  (179) 

So  a  source  doublet  is  the  limit  as  a  source  and  sink  (source  with  negative 
strength)  are  brought  together  with  strengths  inversely  proportional  to  the  dis¬ 
tance  between  them.  The  same  result  is  obtained  by  differentiation  in  equation 
(171). 

In  the  same  manner  in  which  a  source  sheet  was  constructed  by  placing  a  point 
source  in  each  differential  area  of  a  sheet,  we  can  construct  a  doublet  sheet  by 
placing  a  point  doublet  solution  in  each  differential  area  with  the  above  doublet 
solution.  However,  there  were  restrictions  placed  on  the  above  point  doublet  for¬ 
mula.  It  is  limited  to  steady  incompressible  flow.  In  this  report,  doublets  will  be 
formulated  for  unsteady  compressible  flow. 

For  a  point  source,  the  potential  which  arises  at  any  other  point  is  proportional  to 
1/r.  Therefore,  die  potential  field  for  a  source  solution  is  symmetric  with 
respect  to  any  plane  passing  through  the  point  r  »  0,  including  the  (a,  y) 
plane.  It  follows  that  the  potential  field  arising  from  a  source  sheet  in  the  ( x ,  y ) 
plane  is  symmetric  with  respect  to  the  (a,  y)  plane.  If  the  potential  immediately 
above  and  below  the  (a,  y)  plane  is  identical,  then  it  follows  from  equation  (52) 
that  the  pressure  immediately  above  and  below  the  wing  will  be  identical.  There¬ 
fore,  a  single  source  sheet  cannot  generate  a  pressure  difference. 

For  a  point  doublet,  the  potential  is  proportional  to  -z/r3  and  it  follows  that  the 
potential  field  for  a  doublet  sheet  is  antisymmetric  with  respect  to  the  (a,  >•) 
plane.  Since  the  potential  is  antisymmetric,  we  know  from  equation  (52)  that  the 
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pressure  above  and  below  the  (x,  y)  plane  must  also  be  antisymmetric.  So,  in 
contrast  to  the  source  sheet,  we  can  develop  a  pressure  difference  with  a  single 
doublet  sheet  This  is  a  fundamentally  simple  but  important  concept 


SECTION  XI 


The  Acceleration  Potential 


Up  to  this  point,  we  have  been  using  the  velocity  potential  <f>  as  the  unknown 
variable  in  our  linear  aerodynamic  system1.  The  input  to  our  linear  system  is  the 
wing  deformation  h.  The  output  of  our  linear  system  is  pressure  p.  However,  our 
solution  technique  solves  for  the  velocity  potential  first  and  pressure  is  computed 
later  in  a  second  step.  In  this  section,  we  avoid  the  intermediate  step  of  solving 
for  the  velocity  potential  and  solve  for  pressure  directly  with  the  introduction  of 
the  pressure  potential  and  subsequently  the  acceleration  potential2  y. 

We  start  with  the  aerodynamic  potential  equation  (42). 
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Next,  we  differentiate  with  respect  to  t  and  then  x  to  form  the  following  equa¬ 
tions 
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We  multiply  equation  (181)  by  p0  and  equation  (182)  by  pQU  and  add  them 
together  to  obtain  equation  (183). 


1.  St*  William,  Guderiey  and  Lee  for  the  ooe- angular  formulaiioo  for  Sie  potential  which  triaea  from  •  doublet  tbetL 

2.  y  may  elao  be  deacribed  u  •  preture  doublet.  To  be  eontiaiwi  with  ^quadoo  ( 189), «  abould  uea  the  rysbd  p4  to  <Je* 

•crib*  a  pruiuie  doublet.  Htwevts.  we  rwaaitt  coatiiftt&l  with  Vivian  aod  Aadnw*  tod  uae  y .  '* 
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P  Po  (t +  yfx)  „ +  P0  ■ +  „  +  Pc  (♦,  +  UV  „  - 


22 


2  U 
/»2 


Pn<^+^),r 


1 


(183) 


We  recognize  the  following  form  of  equation  (52)  within  equation  (183). 


J»--P,  (184) 

After  making  the  substitution  in  equation  (183),  we  arrive  at  the  pressure  poten¬ 
tial  equation 


ftPxx  +  Pyy  +Pn  “ 


(185) 


Now  the  variable  p  seems  to  be  doubly  defined  as  both  the  pressure  and  the  pres¬ 
sure  potential.  They  mean  the  same  thing.  The  form  of  equation  (185)  is  mathe¬ 
matically  identical  to  the  form  of  the  aerodynamic  potential  equation  (180). 
Only  the  physical  interpretation  is  different  Therefore,  the  elementary  solutions 
to  the  aerodynamic  potential  equation  (180)  are  also  elementary  solutions  to  the 
pressure  potential  equation  (185).  The  elementary  pressure  source  equation  fol¬ 
lows  from  the  elementary  potential  source  equation  (152). 


“  [JKt) 

(186) 

(187) 

of 

[(*-4)2  +  p2(y-q)%pV- 
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(188) 


The  Acceleration  Potential 


We  now  restrict  ourselves  to  harmonics  in  time.  In  other  words,  variable  time 
dependency  is  replaced  by  a  dependence  on  a  constant  frequency.  By  restricting 
the  problem  to  constant  frequencies  co,  we  obtain  great  computational  savings. 
The  computational  cost  of  solving  the  aerodynamic  flow  over  a  body  or  wing  in 
the  time  domain  is  great  in  comparison.  So  equation  (186)  for  a  pressure  source 
ps  with  strength  A  takes  the  complex  form  pseia‘. 

exp(itiit)  =  psemt  (189) 

Now  we  use  the  formula  (187)  for  retarded  time  x  in  equation  (189).  For  the 
modulus  ps  on  the  right  side  of  equation  (189) 


Ps  = 


R 


exp(ia)(t-x))  = 


-  ,  s  A  -/cot  A 
ps(x,y,z)  =  =  -=exp 


1(0 


ur 


(Af  (J.--4)  -R) 


(190) 


Now  ps  is  a  symmetric  function  with  respect  to  the  z  =  t,  plane.  This  means  that 
we  cannot  use  equation  (190)  to  model  a  pressure  difference  across  a  planar 
wing.  We  seek  a  formula  for  a  pressure  doublet.  As  with  the  source  doublet,  we 
differentiate  equation  (189)  with  respect  to  z  to  obtain  the  definition  for  the  pres¬ 
sure  doublet  which  will  be  called  the  acceleration  potential  y  here.  We  divided 
by  pt(  in  order  to  simplify  the  subsequent  formulae  and  to  bring  these  formulae 
in  line  with  the  original  derivations1. 


(  1 1  d 

r  i n 

d  r  -  ,  iau 

'  1  * 

lp.Je  w  = 

[Po. 

.Pa. 

/aw 


f(^,:)e'“"  (191) 

where  the  modulus  of  y  is  the  differentiation  of  equation  (190)  with  respect  to  i. 


V(*,y,r)  = 


ICO 


I  [p*]  = 


(192) 


l.  See  L  V.  Andre vi  end  H  T.  Vjv»sn. 
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Now  \}/  =  is  an  elementary  solution  to  equation  (185)  and  p  =  is  the 
pressure  which  arises  with  the  acceleration  potential.  Some  useful  relations  are 
now  derived.  We  define  the  non-dimensional  pressure  coefficient  Cv  and  imme¬ 
diately  specialize  it  to  the  acceleration  potential. 


r  ..  2P  -  2V 

p  o  U2  IT 
>  0 


(193) 


Next,  we  investigate  the  relationship  between  the  velocity  potential  and  the 
acceleration  potential.  This  is  important  in  deriving  the  boundary  conditions  for 
the  pressure  potential  boundary  value  problem.  We  denote  the  harmonically 
oscillating  potential  as 


few 


<}>  (x,  y ,  z,  /)  =  <j)  (x,  y,  z)  e‘w ‘  (194) 

The  overvar  indicates  the  complex  modulus.  The  overbar  here  does  not  indicate 
the  steady  state  condition  as  used  in  Section  H.  From  \|/  =  —  and  equation 
(184),  we  obtain1 


y(*,y,  z>  0  =  - 

Carrying  out  the  operations 


Wt  +  uTx 


d  1  .-r.  .  iCM. 

(<} )(x,y,z)e  ) 


(195) 


V  (*.  y>  '*)  -  ~U(bx  ( x,  y ,  z)  / co<|>  (x,  y,  z) 


(196) 


We  arrive  at  die  inverse  relation  to  equation  (196)  by  using  p 
equation  (73). 


<)>  (x,  y,  z) 


-1 

* 

r 

"  icoX" 

Uexp 

L  u  J 

j exp 

L  u  J 

— OO 


w  (X^y,  z)dX 


(197) 


1.  Here,  we  can  see  how  the  term  “acceleration  potential”  arises.  It  is  the  total  derivative  of  the  velocity  potential.  The  total 

derivative  is  a  derivative  with  respect  to  time  relative  to  a  steady  translating  frame  of  reference.  The  oruer  of  differentiation  ia 
optional.  Therefore,  the  acceleration  field  is  related  to  the  acceleration  potential  by  the  gradient. 
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The  Integral  Formula 


In  the  previous  section,  we  identified  the  elementary  solutions  to  equation  (185). 
We  derived  the  following  elementary  pressure  doublet  expression  (192),  other¬ 
wise  known  as  the  acceleration  A  ential. 


\j/(j:,y,z)  =  (A)^ 


rxp 


1(0 


(198) 


where  A  =  A  (co)  is  the  amplitude  of  the  oscillations  at  any  given  frequency. 
Also,  we  have 


r  2  2nl/2 

R  =  I  (x-5)2  +  P2(>->1)  +P2(-’-0  J  (199) 

Equation  (198)  gives  the  pressure  field  which  arises  from  a  single  pressure 
doublet  In  this  section,  we  extend  the  single  pressure  doublet  to  a  doublet  sheet. 
As  indicated  in  Section  X  for  the  source  doublet  (velocity  potential),  the 
pressure  doublet  sheet  is  also  suitable  for  modelling  the  pressure  difference 
between  the  upper  and  lower  surfaces  of  a  thin  wing.  Our  goal  in  this  section  is 
to  develop  the  integral  equation1  (224)  which  describes  the  upwash  generated  by 
a  pressure  doublet  sheet.  (For  our  aerodynamic  problem  of  flow  over  a  wing,  the 
upwash  w  is  a  known  function.  In  later  sections  of  this  report,  we  will  see  how  to 
carry  out  the  integral  of  equation  (224)  with  unknown  pressure  Ap.)  First  we 


1.  L.  V.  Andrew*  w>d  H.  T.  Vivian,  pp  15-20. 
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develop  the  relationship  between  the  amplitude  A  and  the  pressure  difference 
A P  across  a  pressure  doublet  sheet  This  is  given  in  equation  (215). 

The  first  step  is  to  carry  out  the  derivative  with  respect  to  z  in  equation  (198).  As 
a  preliminary  step  we  evaluate  the  following. 


8K_  P2(*-S) 
dz  R 

(200) 

a  _  -p2(z- o 

(201) 

Now,  we  carry  through  with  the  derivative  in  equation  (198). 


V  (*»?»*)  =  Ap2(z-0 


'  1  1  ■ 

f" 

_2  ~  —3 

exp 

L  R  R  ] 

1(0 


2  (AT (*-$)-*) 


(202) 


This  is  the  formula  for  a  single  pressure  doublet.  This  formula  describes  the 
pressure  \y  at  coordinates  (x,y,z)  due  to  a  pressure  doublet  at  coordinates 
(4»  0  •  We  now  consider  a  continuum  of  doublets  in  the  (£,,  "H,  ^  =  0)  plane. 

Each  differential  area  dS  is  given  a  doublet  strength  of  AdS.  While  the  choice  of 
differential  partitioning  is  really  somewhat  arbitrary,  we  will  indicate  a  differen¬ 
tial  area  as  a  rectangle  of  area  dS  -  d$dr\  for  the  time  being. 


R 


exp 


1(0 


d%dx\  (203) 


The  next  task  is  to  determine  what  value  takes  as  we  approach  the  surface. 
That  is,  we  evaluate 


lim  V  (x,  y,  z)  (204) 

z-*  0 

—3  —2 

The  presence  of  the  R  and  R  terms  in  die  denominator  of  equation  (203) 
makes  the  integrand  singular  to  order  O  and  O  when  z  =  0  and  when  the 


72 


The  Integral  Formula 


coordinates  (x,  y)  lie  within  the  domain  of  S .  This  can  be  treated  in  the  follow¬ 
ing  manner. 

Construct  a  small  circle  with  radius  po  around  the  point  (x,  y,z  =  0)  in  the 
plane  of  the  wing.  First,  we  integrate  over  the  region  outside  of  the  circle  and 
clearly,  in  the  limit  as  z  0,  the  contribution  to  \y  goes  to  zero.  Only  the 
portion  of  the  doublet  surface  within  the  radius  po  contributes  to  y.  Now,  we 
change  coordinates  such  that 

(*-£)  =  pcos{0) 

(y-‘ n)  =  psin  (0) 


d\dr\  = 


1 

P 


pdpdB  * 


For  p  and  z  sufficiently  small,  we  can  assume 


exp 


) 


»  lim  ( cose  +  /sine)  «  1 
e-*0 


(205) 


__3  „2 

Furthermore,  we  make  z/R  significantly  bigger  than  z/R  by  choosing 
po  and  z  to  be  small.  Finally,  we  assume  A  (£,  T))  is  constant  within  our  small 
circle.  So  we  evaluate 


2rP0 

lim  W  (x,  y ,;)  «  lim  A  |  f 
r-40  .--*o  J  J 

o  o 


(206) 


lim  w  (x,  y,  z) 

X  -4  0 


Iim/4  f 
J 


ol(PV  +  pz) 


3/2 


dp 


(207) 
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lim  \\f  (x,  y,  z)  -  lim  (2jcA)  | 
r-40  i-»0  J 


-P 


r  n  i 

2\3/2 

pyfi+ 

p 

LpJ 

) 

dp 


(208) 


We  now  employ  the  substitution  p  =  pzp. 


limy(x,y,  z)  =  lim  (2tcA)  | 
z->0  ?  0  J 


o  L 


“P 

(l+p2)3/2 


dp 


(209) 


lim  y  (jc,  y,  z)  =  lim  (2tcA)  (p2  +  1) 
i->  o  j  ->  o 


-1/2 


Po  ~  0 


(210) 


We  end  with  the  simple  result  that  as  one  approaches  the  doublet  sheet  from  the 
top  side  that  the  acceleration  potential 

lim  v  (x,  y,  z)  ~  -2nA  (211) 

!-*0 

We  get  a  similar  result  for  the  case  where  the  doublet  sheet  is  approached  from 
the  opposite  side. 

lim  \j7(jc,y,  z)  =  2il4  (212) 

-1-4  0 

The  jump  in  \j/  across  the  doublet  sheet,  going  from  top  to  bottom  is 

A\jr(x,y)  =  4nA  (213) 

We  recall  that  die  relation  between  the  acceleration  potential  and  pressure  is  pro¬ 
portional  to  the  density,  such  that  p  =  \p/p.  So  we  see  from  equation  (213),  the 
pressure  jump  modulus  across  the  doublet  sheet  is  simply 
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Ap(x,y)  =  4npA  (214) 

Our  convention  is  such  that  positive  A p  results  in  positive  lift  with  negative 
pressure  on  the  positive  z  side  of  the  wing  and  positive  pressure  on  the  negative 
z  side  of  the  wing.  So  we  have  the  following  expression  for  A  in  terms  of  Ap. 


(215) 


We  make  the  substitution  into  equation  (203)  and  obtain  the  following  result 


\\f(x,y,z)  = 


R _ 


exp 


i  co 


(M(x-k)-R) 


d%dV[ 


(216) 


We  desire  an  explicitly  linear  relationship  between  the  pressure  jump  across  a 
doublet  sheet  and  the  linear  boundary  condition,  which  is  the  normal  component 
of  flow.  Our  attention  was  diverted  in  order  to  obtain  equation  (215).  We  will  use 
this  relation  later  in  our  effort  to  obtain  the  final  integral  formula.  So  we  reorient 
our  attention  and  begin  our  development  of  this  formulation  with  equation  (197) 
of  the  last  section.  For  a  single  oscillating  doublet  located  at  x  =  wc  ^ave  8111 
oscillating  potential  of 


0  U  y> 


v  -1 

0  a  jjCXp 


'-« 0>  (A  —  4)  ' 

•  7 

r 

"icoA' 

L  U  J 

J  exP 

_1T_ 

y(Ky,z)d\  (217) 


We  substitute  for  y  with  the  expression  found  in  equation  (198). 

-ICO  (x-  £) 


Hxty,z)  =  ~ exp 
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ioA' 

(d 

i 

J  €xp 

.IT . 

rxp 

U 

10) 

4? 


(Af(X-4)-tf) 


\ 

/ 

dX 


(218) 


We  can  move  the  derivative  with  respect  to  z  from  under  the  integral  and  we  can 
combine  exponential  terms.  Equation  (218)  becomes 
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<t >(x,ytz)  = 


-A  1  a 


u  J 


Tz 


exp 


-ia>(x-$) 

U 


x-\ 

lift 


exp 
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{X  MX  R 


u  *p2  *p2JJ 


)dX 


The  z  component  of  velocity  is  related  to  the  velocity  potential  as 

a  - 

w(x,y,z)  =  ^  (<|>) 

Thus,  we  compute  w  (jc,  y,  z)  from  equation  (219). 
w(x,y,z)  = 


(219) 


(220) 


—A' 

a2 

~-ico(*-§) ' 

r  ( 

.  u _ 

dz2 

exp 

L  U  J 

J  rxp 

— oo 

10) 

-  V 

(X  MX  R  \ 

c+^'^2Jj 


dX 


(221) 


We  now  consider  a  continuous  sheet  of  doublets  in  the  identical  sense  as  was 
introduced  in  equation  (203).  Then  substitute  for  A  using  equation  (215). 


(222) 


This  expression  can  be  condensed  slightly1.  Again,  equation  (223)  is  the  formula 
for  the  upwash  w  generated  by  a  harmonically  oscillating  doublet  sheet  in  the 
2  =  0  plane. 


1.  Thereby  pulling  wjunioo  (222)  in  the  identical  fona  uiaibe  rei'erew*  of  Andrew  and  \%vun. 
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w(x,y,z)  = 


-1  1 


4np£/ 


II (A^ 

-1  Cl 


exp 


-ia>  (x-^) 
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a2 

dz2 


rjr-4 


f  i 
J  a 


I© 


(X-MR) 


dX 


d^dr\ 


(223) 


We  choose  to  abbreviate  this  equation  as  follows. 


w(x,y,z)  =  4^1/  JjAp£((x-£),  (y-r\),z)dt>dy\  (224) 


where 


K  (x0,  y0,  z0) 


=  exp( 


-mx0 

U 


(X-MR) 


(225) 


and  we  essentially  repeat  equation  (199). 

*  =  (X2  -t*  p2y2  +  p2z2) 1/2 


(226) 


K  (x0,  y0,  z0)  is  known  as  the  Kernel  function  and  is  the  topic  of  die  next  sec¬ 
tion. 

For  the  purpose  of  this  report,  we  have  restricted  ourselves  to  a  single  planar 
wing.  As  such,  we  may  be  tempted  to  immediately  take  the  limit  as  z  -*  0  in 
equation  (224).  However,  we  still  need  to  evaluate  the  derivatives  with  respect  to 
z  in  equation  (225). 
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The  Kernel  Function 


In  the  last  section,  we  derived  the  integral  equation  (224)  which  relates  the  mod¬ 
ulus  of  die  unknown  pressure  difference  A p  (x,  y )  to  the  modulus  of  the  z  com¬ 
ponent  of  flow  iv  (x,  y,  z) .  We  know  the  z  component  of  flow  at  the  wing 
surface  (boundary  condition)  where  z  =  0.  This  will  manifest  singular  behavior 
which  will  be  addressed  in  the  following  section.  In  this  section,  we  evaluate  the 
Kemal  function. 


K  (x<>>  y<>>  zo)  =  exP 


'  -lOWf.  " 

'  *0 

f 1  r 

exp 

L  u  J 

az2 

ii^[ 

l  CO 


C X-MR ) 


dX 


(227) 


The  objective  of  this  section  is  to  cany  out  the  differentiation  with  respect  to  z  in 
equation  (227).  The  resulting  formula  is  given  in  equation  (257).  We  isolate  the 
integral  expression  in  equation  (227)  and  label  it  as  l0. 


K(*o>y0’  zo > 


exp 


-hox„ 


U  Jaz 


[  U 


(228) 


*• 


1  %exp 

J  R 

«4* 

J 

dX 


(229) 


We  define  two  new  variables  and  Jtl  which  will  be  used  throughout  the 
remainder  of  this  section. 
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2  ,  2v 


'i  =  w+o 


(230) 


*,  = 


This  results  in  the  following  formula  for  R 


(231) 


*  =  (X2+pVJ) 

In  the  evaluation  of  equation  (229),  we  use  a  variable  substitution.  Let 


(232) 


where 


<u  ]** 


Now  use  another  variable  substitution.  We  define  a  new  variable  u. 


u  =  (1  4- v2)  ] 


We  obtain  horn  equation  (234) 


(233) 


4  =  f  ■~T-l--inexp\  (i)  (v^ -Mkx  (v2+  1)  n)  1  dv  (234) 

i.(v  +1)  L  p  j 


(235) 


(236) 


f  exp{~ikxu) 

O  -  j  2  1/2  du 

La+«)  J 


where 


(237) 
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-xd+MR 


(238) 


Equation  (237)  is  used  in  conjunction  with  equation  (228).  In  the  course  of  eval¬ 
uating  the  derivative  of  equation  (228),  we  use  the  chain  rule 


a=ar,j=i3 

<jz 

and  we  obtain 


(239) 
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exp( 
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(240) 


This  can  be  put  in  a  convenient  form  proposed  by  Landahl1.  We  will  eventually 
restrict  our  formula  to  modelling  the  flow  over  planar  surfaces  in  the  (x,  y ) 
plane.  Landahl ’s  equation  represents  a  more  general  case  for  any  planar  surface 
parallel  to  the  x  axis.  We  follow  his  development  and  then  let  z0  go  to  zero. 
Equation  (240)  is  equivalent  to  the  following  equation  (241). 
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Again,  the  expression  for  I0  is  given  in  equation  (237).  Now  we  use  Landahl ’s 
approach  and  make  the  following  substitution  of  variable  in  equation  (237). 

t  »  url  (242) 


=  J 


exp 
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<'i+o 


1/2 


dt 


(243) 


i.  SaatheuticlebyLudahi 
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T  "■  11 


where 


4  =  “44  =  A  (**-*») 


(244) 


In  order  to  evaluate  the  derivatives  of  I0  in  equation  (241),  we  use  Leibnitz's 
Rule1  and  equation  (243). 
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Consider  equation  (245).  From  equations  (232)  and  (244)  we  have 
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Next,  we  carry  out  the  derivative  in  the  integrand  and  equation  (245)  becomes 
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Next,  change  the  variable  of  integration  bach  from  t  to  u  using  equation  (242). 
Equation  (247)  becomes 
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exp(-ikxu) 
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exp(-ikxux) 
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Next,  we  evaluate  the  rest  of  equation  (241).  Starting  with  equation  (247), 


l.  S«e  tl»  fortao*  on  p*p  14  of  thu  ttxi 
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(249) 


Due  to  the  complexity  of  the  formulation,  we  partition  the  right  hand  side  of 
equation  (249). 
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Hx  and  H2  are  defined  in  equations  (251)  and  (254). 


Now,  use  Leibnitz’s  rule  in  equation  (251). 
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('I  +  O 


f  —IQ)/  "j 

1  \ 

/ 

77 

exp 

,  — - r.^ 

5/2 


•icori 


\\ 


V  J 

U»t  +  4>S/aJ 


*i 

^ri 


(252) 


We  get  the  following  result  using  die  relationship  between  t  to  u  in  equation 
(242). 


expi-ikiU) 


exp(-iktu)  ^ 

,  2.3/2 
(1  +ii)  ; 


Next,  we  evaluate  //2  in  equation  (250). 


(253) 
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H2~&1 


-M 

L  R  . 


exp 


-ia >t 


U 


2S  i/2 


(254) 


We  cany  out  the  differentiation  operation. 


i(Q~ 

r*.i 

M- 

—/cor.  * 
e*P(  v  ) 

.77. 

3* I . 

. R . 

1 - 

/-s 

+ 

mV* 

bJ 

\ _ 

(255) 


Now  use  equation  (246)  in  equation  (255).  After  some  basic  algebraic  operations 
the  following  result  is  obtained. 


H , 


f'MfY 

(-i*,#,) ' 

l  rxR2  ) 

U’JJ 

j  *  ,  2.1/2 

L  (1+«i)  J 

(256) 


WSJ  V r,R-  J 


expHkjU,)  ~j 
/i  2.  3/2  I 

(i+«p  J 
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We  now  have  an  algebraic  expression  for  the  Kernel  function  given  in  equation 
(241)  using  equations  (248)  and  (250).  Equations  (253)  and  (256)  are  used  in 
place  of  equation  (250).  It  is  left  to  the  reader  to  consolidate  these  equations.  The 
operation  is  simple;  it  is  not  practical  to  display  such  a  lengthy  formula  in  this 
text  This  completes  our  differentiation  with  respect  to  z  in  equation  (727). 


Landahl  has  provided  a  compact  form  of  our  lengthy  formula.  At  this  point,  we 
switch  over  to  Landahl’s1  notation  in  order  to  follow  his  work  closely.  We  will 
use  the  expanded  terms  following  equation  (241).  Equation  (241)  can  be  placed 
in  the  form 


[K1Tl+K2T2] 

K(x0,y0,z0)  =  exp  - 2 - 


(257) 


where  the  terms  are  directly  relatable  to  the  two  terms  of  equation  (241). 


(258) 


(259) 


In  our  development,  we  have  limited  ourselves  to  a  doublet  sheet  in  the  (x,  y) 
plane.  The  formula  f  iven  by  Landahl  is  more  general,  representing  any  doublet 
sheet  which  is  parallel  to  the  x  axis.  We  make  die  appropriate  modification  tc  the 
expressions  for  Tx  and  T2  by  setting  y(s)  =0  and  y(a)  =  0.  Thus  T1  and  T2 
become 


(260) 


1.  Tike  note  tbat  in  the  report  by  Landahl,  there  is  »n  error  in  the  K2  term.  Thii  u  clear  when  compared  to  equation  (241)  in 
this  section. 
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Ti 


(261) 


We  use  equation  (248)  in  Landahl’s  expression  for  Kx  and  we  use  equation  (250) 
in  his  expression  for  Kv  We  obtain  Landahl’s  resulting1  formula  for  Kx  and  K2. 


~Mrx~ 

exp  (-ikxux) 

L  R  J 

M  2,1/2 

L  (1  +  wi)  J 

(262) 


3  /2  + 
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e-H  : 
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exp  (~ikxux)  1 

L  R2  J 
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N 

+ 

R 

>—  K> 
w 

t-A 
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+ 


+  (2)  + 


'Mrxux~ 
.  R  . 


exp(-ikxux) 
/ 1  ,  2,  3/2 


(263) 


where 


and 


f  (exp  (~ikxu)  \ 


J  t 

u.  N 


M  .  2,3/2 

(l+«  ) 


du 


(264) 


'2 


/exp  (~ikxu)  \ 

,(l+u2)5/2J  “ 


(265) 


Differentiation  with  respect  to  z  within  the  kernel  function  equation  (227)  is 
complete.  The  resulting  equation  for  the  kernel  function  is  given  by  equation 


1.  The  diflerenoe  in  the  sign  between  this  reeult  and  Lindahl's  remit  is  accounted  for  in  the  minus  sign  added  in  equa¬ 
tion  (224). 
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(257).  We  might  feel  free  to  take  the  limit  as  z  goes  to  zero  for  our  planar  wing. 
In  equation  (257)  it  turns  out  that,  because  of  the  singular  nature  at  z  =  0,  this  is 
not  entirely  recommended.  However,  we  note  that 


lim  =  o 


zo->0 


except  when  y0  =  0  in  which  case 


lim  T 
*0  ->  O' 1  2 


=  1. 


We  also  note  that  K2  is  finite  everywhere.  When  K2T2  is  added  to  KlTl  in  equa¬ 
tion  (257)  and  K  is  integrated  in  equation  (224),  the  contribution  of  K2T2  is  zero 
when  z0  is  zero.  Therefore,  for  a  planar  wing,  we  can  immediately  state, 


KmK(x0,y0,z0) 


<*o 

'-i(0xo- 

exp 

L  u  J 

(266) 


Equation  (264)  for  Ix ,  can  be  modified  for  improved  computation.  We  assume1 

ux  ;>o. 


r  _  [exP(-ik\u) 
7i  ~  J  7372 
«i  (l  +  «  ) 

Integrate  by  parts  to  obtain 


du  »  J  [exp  (-ik^)] 


du 


L(i +« ) 


3/2 


A- 


XP (-'*!“)  { - i7T7i] 

Ml  +  « )  >J 


ik\  ( - “ 

uWl  +  U 


1/2 


exp  {~ikxu)  du 


^(1  +  0  ' 


(267) 


(268) 


1.  Use  equation  (275)  for  <  0. 


87 


T^^pjeJjFunctioo 


Evaluate  the  term  in  square  brackets. 


h  = 


exp  (-ikx°°)  -  exp  (— i&jiij) 


Mi 


(1  +  u]) 


1/2 


*/(— 


1/2 


exp(-ikxu)du 


(269) 


The  term  exp  (-Z&joo)  is  somewhat  meaningless.  We  fold  this  back  in  the  inte¬ 
gral 


Ix  =  -exp(-ikxu  j) 


(1 


,  1/2 

+  u])  ) 


ik\  (-  1  + 


M  ^ 


1/2 


exp{-ikxu)du  (270) 


«,  v  ( 1  +  u  )  7 

This  can  be  abbreviated  as  follows. 


h  =  exp(-ik !«!) 


1  -  f  U1 

+  HVi) 

L  Ml+nl)  J 

(271) 


where 


=  exp(ikxux)  J 


“iL 


1- 
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(1+M2) 


1/2 


exp  (- ikxu )  du  (272) 
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Laschka1  provides  the  following  accurate  approximation  for  u  £  0 

“  X  a"e~°C“  (273) 

n  =  1 

where  c-0372  and  the  an  are  given  in  the  following  table 

n  an 

1  +0.24186198 

2  -2.7918027 

3  +24.991079 

4  -111.59196 

5  +271.43549 

6  -305.75288 

7  -41.183630 

8  +545.98537 

9  -644.78155 

10  +328.72755 

11  -64.279511 

Substitute  equation  (273)  into  equation  (272)  to  obtain  an  approximate  expres¬ 
sion  for  Jx 


1- 


u 


\ 


1/2 


■d+o  n 


H 


Jx « 


n=  lL 


~ncu 


2  2  .  ,2  . 

n  c  +  kx  J 


(nc~ikx) 


(274) 


We  use  equation  (274)  in  equation  (271)  to  obtain  a  formula  for  Ix .  This  formula 
is  valid  only  for  ux  >  0.  We  see  that  the  integrand  of  Ix  in  equation  (264)  is  sym- 


1.  S<e  L«»chk«  u  pointed  out  by  Glaring  et  aL  on  page  55  of  Part  I,  VjI  L 
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metric.  For  Wj  <0  we  can  take  advantage  of  the  symmetry  of  the  integrand  and 
still  use  the  algorithm  of  equation  (274).  We  evaluate  the  real  and  imaginaiy 
components  separately.  For  ul  <  0 : 

h  ( uv  kt)  =  2 Re  [7j  (0,  *j)  3  -  Re  [Ix  (-»„  kt)]  +  Um  [Ix  (-uv  kx)  ]  (275) 
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The  Doublet  Lattice  Method 


There  is  no  precise  definition  of  the  doublet  lattice  method  and  the  associated 
formulae.  Basically,  this  Section  restricts  the  approach1  of  Giesing  et  al.  to 
planar  wings.  In  addition,  we  deviate  in  the  treatment  of  the  sweep  angle.  This 
deviation  will  be  identified. 

For  a  wing  in  the  z  =  0  plane,  we  have  the  following  integral  formula  from 
equation  (224). 

w  (x,  y ,  0)  =  ((*“§)>  (y  -  “H) » 0)  d%dr\  (276) 

with  the  following  supplementary  formulae  from  equations  (266),  (262),  (264), 
(231)  and  (238).  Here,  we  have  substituted  €  for  z  to  emphasize  the  limitation 
process. 


K(xo>yo>0)  =  hm 


(277) 


M\y0 

exp(-ikiUx) 

L(4+P2^)  J 

,<  ,  .  2% 1/2 

(278) 


1.  See  Gietiog,  Ktlmtn  and  Rodden. 
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h 


exp  (-i^u) 
(1  +  u2)3/2 


du 


(279) 


(280) 


u 


l 


M{x20^2yl)m-x0 

wp2 


(281) 


—2  2 

The  Kernel  function  in  equation  (277)  is  0  singular  as  y0  goes  to  zero.  This 
requires  special  consideration  as  one  integrates  equation  (276).  This  singularity 
occurs  when  the  y  and  T|  coordinates  are  the  same.  Furthermore  Kx  is  singular 
when  x0  and  y0  are  both  zero.  This  occurs  when  (x,  y)  is  coincident  with 
(£,  q) ,  Finally,  as  ux  ranges  from  -<*  to  +<*>,  we  must  give  special  consider¬ 
ation  to  equation  (279).  The  following  question  remains.  Can  we  evaluate  the 
integral  equation  (276)  with  these  mathematical  difficulties.  All  these  difficulties 
can  be  overcome  analytically  if  one  uses  approximating  functions.  In  this  ana¬ 
lytic  approach,  one  appeals  to  principle  values1.  The  principle  values  procedure 
requires  the  limit  as  e  -» 0  be  taken  as  the  absolutely  final  step.  If  one  is  ever  in 
doubt  about  this  procedure,  it  may  prove  reassuring  to  evaluate  the  integrand  for 
several  cases  and  plot  the  value  as  the  integrand  approaches  the  limiting  singu¬ 
larity. 

The  doublet  lattice  method  is  an  empirical  device  which  simplifies  the  integra¬ 
tion  of  this  singularity.  But  primarily,  the  advantage  gained  by  the  doublet  lattice 
method  is  the  relative  simplicity  in  the  resulting  computer  program,  especially 
for  complex  configurations.  There  are  other  methods  which  are  not  as  simple  to 
implement  On  the  other  hand,  there  is  a  simpler  method  called  the  doublet  point 
method2.  With  the  doublet  lattice  method,  the  continuous  pressure  doublet  sheet 


l.  See  Mangier 
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in  equation  (276)  is  replaced  by  a  set  of  pressure  doublet  lines  with  finite  length. 
In  figure  3  we  picture  the  pattern  of  nine  doublet  lines  for  a  rectangular  wing. 
Each  line  is  contained  in  its  own  box.  The  simplicity  is  that  all  boxes  are  treated 
identically,  regardless  of  its  proximity  to  the  wing  boundary  (i.e.  leading  edge, 
trailing  edge  or  wing  tip).  Other  methods,  based  on  a  continuous  doublet  distri¬ 
bution  require  special  square  root  singularities  in  the  pressure  distribution  near 
the  wing  boundaries. 


Figure  3.  A  Rectangular  Lattice 

The  doublet  line  is  placed  at  the  quarter  chord  of  each  box.  (To  call  this  a  “dou¬ 
blet  lattice”  is  really  a  misnomer.  If  one  views  the  doublet  line  segments  alone, 
no  lattice  is  formed.  The  name  “doublet  lattice"  arises  from  the  conectly  named 
“vortex  lattice”  methods1  applicable  to  unsteady  incompressible  or  steady  com¬ 
pressible  flow  over  planar  wings.)  The  upwash  w  (*,  y,  0)  is  evaluated  at  the  3/4 
chord  midspan  of  each  box.  The  empirical  nature  of  the  doublet  lattice  method 


2.  See  Ueda  and  Dowell 
1.  See  June* 
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arises  in  the  choice  of  the  3/4  chord  and  1/4  chord.  By  no  means  is  there  any 
mathematical  proof  that  this  is  the  correct  location.  As  a  matter  of  fact,  for  non 
rectangular  wings  with  swept  and  tapered  boxes,  the  programmer  is  left  to  his 
own  devices1  to  invent  meaningful  1/4  chords  and  3/4  chords.  However,  in 
defense  of  the  doublet  lattice  method,  this  level  of  empiricism  is  probably  con¬ 
sistent  with  the  level  of  approximation  employed  in  formulating  the  linear  aero¬ 
dynamic  potential  equation  (42).  One  must  be  on  guard  and  realize  that  the  flow 
field  generated  by  this  lattice  of  doublets  will  not  be  smooth,  especially  near  the 
wing  surface.  What  is  important  is  that  the  upwash  at  the  3/4  chord  is  approxi¬ 
mately  the  same  whether  one  has  a  constant  strength  doublet  line  at  the  1/4  chord 
or  has  a  continuous  doublet  sheet  with  the  correct  strength. 

Rather  than  dwell  on  rectangular  wings,  we  assume  we  can  invent  a  meaningful 
location  for  the  doublet  line  and  the  upwash  point  for  general  trapezoidal  boxes. 
A  discretized  swept  and  tapered  wing  is  pictured  here  in  Figure  4. 


Figure  4.  A  Swept  and  Tapered  Lattice 


1.  For  in«l*nce,  Gieaing  et  al.  u»ed  *  hybrid  approach.  In  their  calculation*.  They  itait  with  a  swept  doublet  line.  However, 
the  re  Rilling  integral  formula  u  carried  out  over  an  uarwcpt  doublet  line.  Thia  point*  out  the  cmpiricura  of  the  doublet  lattice 
method. 
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So  the  area  integral  of  equation  (276)  is  truncated  to  a  line  integral  along  the  1/4 
chord  of  each  box.  This  line  is  depicted  in  Figure  5.  Furthermore,  we  assume  A p 
is  spatially  constant  for  each  box.  For  a  rectangular  box,  the  box  chord  is 
denoted  as  A£.  Equation  (276)  becomes 


w (x, y, 0)  =  ^^-limjA' ((x-Z,),  (y-r\),e) dl  (282) 


4np  U  e 


Figure  5.  Local  Swept  Coordinates 


Substituting  for  K  from  equation  (277)  gives 
w  (x,  y,  0)  = 


-ApA^,.  <y — ) ) e-tp [  ——jj — 


AnpU  Je 


(y-‘ n)2+e2 


(283) 
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Let  /j  — L  and  l2  -  L.  We  abbreviate  equation  (283)  as  follows 


w  (x,  y,  o) 


--ApA^ 

lim  f 

.  4 npU 

e->0J 

-L 

(y-‘ n)2+e2 

dl 


(284) 


where 


K  (xQ,  y0, 0) 


Ki  (xo>y0)  exP 


--mx0- 
.  U  . 


(285) 


Now,  it  is  clear  K  is  a  complex  function.  It  turns  out  that  K  can  be  approximated 
with  a  complex  parabolic  function  of  /. 


K(x0,y0)  -  A0-^All^A2t1 

where  A0,  Al  and  A2  are  complex  coefficients.  We  identify  the  coordinate 
(xL>yi)  *°  represent  (x0,y0)  at  /  ==  —L.  At  the  opposite  end,  we  identify  the 
coordinate  (xRtyR)  to  represent  (x 0>y0)  at  /  a  L.  Finally,  at  the  midpoint,  we 

use  (*o?c)  t0  represent  (x0,yo)  at  /  =  0.  Equation  (286)  can  be  formulated 
as 


R  (x0,  y0) 


Hl-L) 
.  1L2 


&  (xL,  yL) 


K(xc,yc)  + 


r/q+0 

.  2L2 


%(xR>yn) 


(287) 


which  is  easily  verified  at  /  =  -L,  /  =  0  and  /  =  L.  We  regroup  equation  (287) 
in  terms  of  common  powers  of  /.  Thus,  we  can  identify  the  coefficients  in  equa¬ 
tion  (286)  as  follows 
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Aq  -  K(xc,yc ) 


(288) 


= 


k  (xr>  yR)  k  (xl>  yD 
_ 


(289) 


^2  “ 


K  (xD  yL)  -  IK  (x0  yc)  +  K  (xR,  yR) 
2L 2 


(290) 


Now  we  substitute  equation  (286)  into  equation  (284).  Furthermore,  we  note 
rj  =  /sin  A. 


w’Uy,  0)  = 


'  -ApA^ ' 

L 

lim  f 

A0+All+A2tl 

_  4np(/ 

6  -40  J 
-L 

_  0?~/sinA)2  +  e2. 

dl 


(291) 


We  abbreviate  equation  (291)  in  the  following  fashion. 


wtx.y.O)  at 


-ApA§ 

JnpU 


{Bq+Bi  +  Bz) 


(292) 


The  definitions  of  B0 ,  Bx  and  B2  are  given  as  equations  (293),  (297)  and  (299). 


Bq  -  Urn  [ 
U  £->0J 


L  (sin A )2/2-  (2>sinA)/-f  (/  +  e2)  j 


1  dl 


(293) 


We  integrate  equation  (293)  to  obtain  the  following  inverse  tangent  function 


Bn  -  lim 

u  5— >  o 


esinA 


a  tan 


/sinA— yl  L 


(294) 
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We  use  the  standard  tangent  identity  for  the  difference  of  two  angles  to  simplify 
equation  (294).  We  obtain 


JiSoLesinAj 


atan 


2eLsinA 


|_e2  +  y2-L2(sinA)2 


Now  we  take  the  limit 


(295) 


*0  = 


2 LA, 


y1-L1(  sinA)2 


(296) 


We  follow  a  similar  procedure  for  to  obtain  algebraic  expression  for  Bx  and  B2- 
The  definition  of  Bx  is 


j3,  =  lira  f 
1  e-+0J 


AJ 


U  (sinA)2/2 -  (2ysinA) /  +  (y2  +  e2) 


d\ 


After  integrating,  we  take  the  limit  and  obtain  the  following  formula. 


B ,  = 


r  ] 

Inn 

(sinA)2/!,2-  (2ysinA)£  +  (y2) 

.2  (sinA)2. 

(sinA)2L2+  (2ysmA)£.4-  (y2) 

yA, 

r  2  l  i 

sinA 

y2-L2(  sinA)2. 

The  definition  of  £•,  is 


(297) 


(298) 


B 2  = 


L 

lim  f 
t-»0j 
-L 


AS 


dl 


(sinA)2/2 -  (2ysinA)/  +  (y2  +  e2) 

Again,  after  integrating,  we  take  the  limit  and  obtain  the  following  formula. 


(299) 
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2 LA2  ' 

+ 

'2  yA2 

Inn 

(sinA)2L2-  (2ysinA)L  +  (y2) 

( sin  A) 2 

sinA_ 

lug 

_  (sinA)2L2+  (2ysinA)L  +  (y2)  _ 

y\ 

r  2  L  I 

( sin  A) 2 

y2  -  L2  ( sinA) 2  _ 

(300) 


The  values  for  f?0,  Bx  and  B2  are  substituted  in  equation  (292)  to  obtain  the  rela¬ 
tionship  between  A p  of  one  element  and  the  w  generated  at  the  control  point  of 
another  element  Again,  AS,  is  the  average  chord  of  an  element  box. 


It  is  important  for  the  doublet  lattice  user  to  understand  the  approximations 
incurred  in  discretizing  a  doublet  sheet  into  trapezoidal  boxes.  It  should  be 
immediately  obvious  that,  since  we  have  assumed  the  pressure  to  be  constant 
within  each  box,  a  sufficient  number  of  boxes  is  required  to  capture  the  steady 
state  (zero  frequency)  pressure  function  accurately.  It  is  not  obvious  that  we  need 
to  increase  the  number  of  boxes  as  we  increase  the  frequency  of  oscillation.  For 
instance,  the  pressure  field  over  a  rigid  wing,  plunging  at  high  frequency,  is  not 
trivial  and  requires  a  significant  number  of  boxes  to  resolve  the  standing  (pres¬ 
sure)  waves.  The  required  box  density  depends  on  a  combination  of  wing  defor¬ 
mation  and  die  frequency  of  motion.  The  box  density  should  be  increased  as  the 
deformation  becomes  more  spatially  wavy  and  as  the  temporal  frequency  of 
motion  increases  The  doublet  lattice  user  must  perform  convergence  studies  to 
determine  the  appropriate  box  density  for  their  application. 
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The  purpose  of  this  section  is  to  introduce  a  clear  and  simple  version  of  a  doublet 
lattice  computer  code.  The  mathematics  of  this  text  are  tedious.  Unless  one  is 
somehow  inspired,  these  mathematics  seem  to  exceed  the  bounds  of  reasonable¬ 
ness.  It  is  possible  that  one  may  overcome  this  hurdle  by  browsing  through  a 
well  annotated  version  of  a  doublet  lattice  code.  Other  existing  computer  codes 
for  full  aircraft  configurations  are  difficult  to  follow  because  of  the  programming 
details.  Clearly,  this  is  not  a  criticism  of  the  usefullness  of  these  codes.  Afterall, 
the  aerospace  community  has  depended  on  them  for  over  twenty  years  now. 
They  function  well  for  a  wide  variety  of  configurations. 

At  the  beginning  of  this  effort  to  provide  a  tutorial  on  the  doublet  lattice  method, 
this  author  had  a  hope  that  he  could  start  with  an  existing  code,  simplify  it  to  the 
planar  case  and  then  add  comments.  It  turned  out  that  it  was  more  effort  than  was 
warrented.  This  author  decided  that  the  algorithm  was  so  conceptually  simple 
that  he  would  develop  iris  own  code.  Afterall,  the  whole  raison  d’etre  for  the 
doublet  lattice  method  is  that  it  is  relatively  easy  to  encode  on  a  computer. 

The  choice  of  computer  language  was  not  easily  made.  While  there  is  a  tremen¬ 
dous  sentiment  for  engineers  to  use  FORTRAN,  it  is  nowhere  near  the  most 
overall  popular  computer  language.  The  C  language  is  very  popular,  especially 
on  personal  computers  (PC),  and  it  is  adequate  for  encoding  the  doublet  lattice 
method.  Hie  most  important  feature  of  the  C  language,  as  far  as  this  tutorial  is 
concerned,  is  the  easy  to  read  format  Comments  can  be  placed  just  about  any¬ 
where.  It  is  much  easier  to  point  out  the  relationship  between  the  lines  of  the 
computer  code  to  the  material  of  this  text  In  addition,  the  C  language  is  very 
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popular  within  the  computer  graphics  community.  With  the  doublet  lattice 
method  encoded  in  C,  it  is  far  easier  for  the  PC  programmer  to  connect  it  with  a 
graphics  library.  This  author  believes  that  some  day,  a  visually  enhanced  version 
of  a  doublet  lattice  code  on  a  PC  computer  will  be  used  to  effectively  motivate 
students  toward  the  study  of  unsteady  aerodynamics.  The  main  advantage  of  the 
FORTRAN  language  is  the  COMPLEX  data  type.  There  is  no  equivalent  data 
type  in  C.  This  author  chose  to  accept  this  shortcoming  and  use  the  C  language. 

Tbs  source  code  for  this  doublet  lattice  code  for  a  simple  trapezoidal  wing  is 
given  in  Appendix  A.  The  comments  within  the  source  code  are  sufficient  for 
one  to  relate  to  the  equations  of  this  report  The  example  input  and  example  out¬ 
put  arc  given  in  Appendicies  B  and  C.  This  example  case  is  for  a  simple  plung¬ 
ing  rectangular  wing  of  aspect  ratio  two. 

As  mentioned  earlier,  there  is  no  strict  mathematical  basis  for  the  doublet  lattice 
method.  It  seems  to  work  for  rectangular  wings.  For  swept  and  tapered  wings, 
we  define  what  we  mean  by  the  1/4  chord  and  the  3/4  chord.  Geising  et  al.  seem 
to  use  a  hybrid  approach.  The  integration  process  treats  the  line  doublet  as 
though  it  was  not  swept.  However,  the  Kernel  function  is  evaluated  at  three 
points  along  the  swept  doublet  line.  The  approach  taken  in  Appendix  A  is  to 
integrate  along  a  swept  doublet  line.  The  effect  of  changing  the  sweep  of  the 
doublet  line  may  be  an  interesting  topic  for  study. 

In  a  sense,  the  doublet  lattice  method  is  empirical.  Giesing  et.  al.  point  out  that 
for  steady  state  analysis,  the  integral  formulae  can  be  integrated  “exactly”.  In 
order  to  achieve  increased  mathematical  accuracy  (this  does  not  guarantee  that 
correlation  with  test  data  will  improve)  Giesing  et.  al.  chose  to  subtract  out  the 
steady  state  component  computed  by  the  doublet  lattice  method  and  replace  this 
component  with  “exact”  computations.  While  this  could  be  done  in  Appendix  A, 
it  wasn’t  The  point  of  Appendix  A  was  to  explain  the  implementation  of  the 
doublet  lattice  method  without  added  complication.  The  reader  should  be  able  to 
see  how  to  implement  a  correction  to  the  steady  state  component  However,  this 
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author  is  not  ready  to  say  that  anything  is  really  gained  for  the  effort.  This  is 
another  suggested  topic  for  study. 

All  the  equations  presented  in  this  report  made  no  mention  of  the  units  of  mea¬ 
sure.  It  turns  out  that  there  may  be  a  benefit  in  non-dimensionalizing  these  equa¬ 
tions.  The  non-dimensional  solution  depends  on  the  Mach  number,  the  reduced 
(non-dimensional)  frequency,  the  shape  of  the  wing  planform  and  the  non- 
dimensional  deformation.  The  solution  can  then  be  scaled  to  meet  a  variety  of 
different  conditions  such  as  vehicle  velocity  and  air  density.  Therefore,  at  the 
expense  of  possibly  complicating  the  interpretation  of  Appendix  A,  all  the  vari¬ 
ables  are  assumed  to  be  non-dimensional.  Non-dimensional  time  t  is  scaled  by 


and  non-dimensional  length  x  is  scaled  by  the  characteristic  length  b 


(301) 


X  =  g  (302) 

The  reader  should  have  no  trouble  in  fomulating  non-dimensional  up  wash  w. 
One  merely  divides  the  dimensional  upwash  by  the  frees tream  velocity  U. 

The  example  case  is  for  a  simple  rectangular  wing  with  aspect  ratio  of  two.  Only 
half  of  the  wing  is  modelled.  The  wing  is  symmetric  about  the  x  axis  in  all 
respects.  The  wing  is  plunged  with  a  reduced  frequency  of  one.  There  can  be  no 
correction  for  steady  state  because  the  zero  frequency  load  is  zero.  The  solution 
agrees  with  data  computed  with  the  method  of  Giesing  el  al. 

The  point  of  this  report  is  not  to  provide  a  detailed  explanation  of  the  implemen¬ 
tation.  The  point  is  to  compile  all  the  mathematics  which  lead  to  the  doublet  lat¬ 
tice  method  in  one  single  document  This  has  been  done.  The  utility  of  the  code 
in  Appendix  A  is  not  assured.  The  author  of  this  report  decided  to  include  this 
code  with  the  hope  that  its  mere  inclusion  would  help  illuminate  the  doublet  lat¬ 
tice  method.  The  reader  should  feel  free  to  use  the  example  source  code  pro¬ 
vided  here  as  a  starting  point  for  developing  their  own  utility.  However,  before 
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doing  so,  one  should  give  serious  consideration  to  using  the  code  erf  Geising  et 
al.  As  indicated,  their  code  is  very  versitile  and  is  well  proven.  It  has  been  the 
mainstay  of  aeroelastic  analysis  and  design  for  two  decades. 
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APPENDIX  A 


The  Doublet  Lattice  Program  Source  Code 


/*  the  following  parameters  may  be  adjusted  */ 

♦define  MAX_DIV_X  20/*  maximum  number  of  chordwise  boxes  */ 
♦define  MAX_DIV_Y  20/*  maximum  number  of  spanwise  boxes  */ 
♦define  MflX_POLY  20/*  maximum  number  of  polynomial  coeff  */ 
♦define  PAUSE_ON_OUTPUT  2  /*  time  to  pause  for  reading  output  */ 
/*  Remember  to  recompile  after  adjusting  the  above  */ 

/*  do  not  adjust  the  following  parameters:  */ 

♦define  MAXBOX  (M6X_DIVJ{  *  M&XJ3IVJO 
♦define  MAXDIM  (MAXBOX  *  MAXBOX) 

♦define  ABS(x)  ((<x)<0)  ?  -  (x)  :  (x) ) 

♦define  PI  (3.141592653589793) 

♦define  EPS  (1.0e-6) 

♦define  BIGP  (l.Qe+20) 

♦define  BIGM  <-1.0e+20) 

/*  end  of  define  */ 


struct  element 

( 

float  xi,yi;/*  coord  of  inboard  1/4  chord  */ 

float  xm,ym;/*  coord  of  midspan  of  1/4  chord  */ 

float  xo,yo;/*  coord  of  outboard  1/4  chord  */ 

float  xc,yc;/*  coord  of  control  point  at  3/4  chord  */ 

float  chord, area;/*  box  chord  and  area  */ 

float  xcent, ycent; /*  x  and  y  coord  of  centroid  */ 

In¬ 
struct  trapezoid 
{ 

int  aymm; 

int  nucs_box_x,  nua_boxjy,  total_boxes; 
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float  xible,yible; 
float  xibte,yibte; 
float  xobte,yobte; 
float  xoble,yoble; 
float  mean_chord, area; 
>/ 

struct  polynomial 
{ 

float  a; 
int  px,py; 

} ; 


/★********************************************i******************* 

*  This  is  a  doublet  lattice  code  for  a  single  trapezoidal  wing.  * 

*  You  can  assume  symmetry  or  anti  "-symmetry  about  the  x  axis.  * 

*  This  code  was  written  by  Max  Blair  of  USAF  Wright  Labortory.  * 

*  Neither  Dr  Blair  or  the  USAF  assume  legal  responsibility  for  * 

*  potential  errors  which  exist  in  this  computer  code.  The  user  * 

*  is  encouraged  to  validate  the  code  for  his  or  her  design  cases* 

*  of  interest.  Send  comments  to:  * 


* 

*  Dr  Max  Blair 

*  WL/FIBRC 

*  WPAFB,  OH  45433“- 6553 

* 

*  This  code  was  written  primarily  for  educational  purposes.  For 

*  complete  aircraft  configurations,  the  user  is  encouraged  to 

*  use  the  doublet  lattice  codes  H7WC  and  N5KA. 

* 

*  The  input  is  placed  in  files  dl. INPUT  and  be. INPUT 

* 


*  dl. INPUT: 


*  BLAIRCRAFT 

*  5.0 


* 

0.5 

* 

1.0 

* 

3 

* 

0.0 

0.0 

* 

10.0 

0.0 

* 

10.0 

10.0 

* 

0.0 

10. G 

• 

10 

* 

10 

2100  ATTACK  FIGHTER(title  line) 
characteristic  length  (b) 

Mach 

reduced  frequency  wb/U 

s:  symmetric  a:  anti-symmetric  n:  no  symmetry 
x  and  y  coord  of  inboard  leading  edge 
x  and  y  coord  of  inboard  trailing  edge 
x  and  y  coord  of  inboard  leading  edge 
x  and  y  coord  of  outboard  leading  edge 
number  of  chordwise  cuts  (discretized  x) 
number  of  spanwise  cuts  (discretized  y) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 
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* 


* 


*  be. INPUT: 

*  flag  constant  x  power 

*  1  -1.0  0 

*  0  -1.0  1 

*  0  -1.0  0 

*  0  -1.0  2 

*  0  -1.0  1 

*  0  -1.0  0 

*  end  of  data 


y  power 
0 
0 
1 
0 
1 
2 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


*  interpretation  of  be. INPUT:  * 

*  w(x,y)  «  a00  +  al0*x  +  a01*y  +  a20*xA2  +  all*x*y  +  a02*yA2  * 

*  where  w  has  been  non-dimensionalized  by  the  velocity,  U.  * 

*  Only  lines  with  "1"  in  the  first  column  is  considered  data.  * 

*  Replace  the  "1"  with  a  "0*  to  ignore  any  data.  * 

*  A  line  which  begins  with  an  "ew  will  terminate  the  input.  * 

*  There  must  be  at  least  one  line  which  begins  with  an  M*.  * 

*  * 


*  NON-DIMENSIONAL  INPUT:  * 

*  If  wing  coordinates  are  already  normalized  with  respect  to  * 

*  a  characteristic  length,  then  input  b*l.  * 

*  * 


*  DIMENSIONAL  INPUT:  * 

*  If  wing  coordinates  are  input  in  other  units  (such  as  inches)  * 

*  then  input  any  value  for  b  such  as  the  mean  aerodynamic  chord  * 

*  in  consistent  units  (inches) .  The  upvash  is  input  in  * 

*  non-dimensional  form,  normalized  with  respect  to  the  free  * 

*  stream  velocity.  Non-dimensional  pressure  coefficient  will  * 

*  be  printed  out.  * 

*  * 


*  Output:  * 

*  complex  modulus  of  the  pressure  coefficient  (Cp)  at  each  box  * 

*  Cp  »  pressure/ (density*velocity  squared)  * 

********************** ******************************* ***********/ 


♦  include  Ml. define" 
♦include  <math.h> 
♦include  <atdio.h> 
♦include  *di. structure" 


main  () 

( 

FILE  *fopen  () ; 

FILE  *aicdat,  *odat;/*  pointers  to  I/O  files  */ 
int  discretize  < ) , read^input ( ) , Kbar ( ) , be ( ) ; 
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unsigned  seconds; 

int  i,symm, ic; 

int  ierr;  /*  error  code  */ 

int  rb,sb;  /*  receiving  box  and  sending  box  indicies*/ 

int  vector_index;  /*  how  the  [D]  matrix  is  placed  in  a  vector  */ 

float  xO,yO;  /*  distance  from  receiving  to  sending  coordinates  */ 

float  M,  k;  /*  Mach  and  reduced  frequency  */ 

float  b,  b2;  /*  characteristic  length  (also  bA2)  */ 

float  ul,beta2;  /*  ul  is  eqtn  281  */ 

float  KbarrL,KbariL;  /*  real  and  imag  Kbar  for  left  point  */ 
float  KbarrC, KbariC;  /*  real  and  imag  Kbar  for  center  point  */ 
float  KbarrR, KbariR;  /*  real  and  imag  Kbar  for  right  point  */ 
float  L,L2,dLx,dLy;  /*  related  to  length  of  the  doublet  line  */ 
float  Yt,Yt2;  /*  y  relative  to  the  sending  midpoint  */ 

float  sl,sl2,sl3;  /*  sine  of  the  doublet  line  angle  (fig  £)*/ 
float  factor,  facto,  fact ljr  f 9 ct2,  fact 3;  /*  working  space  */ 
float  AOr,  AOi,  Air,  Ali,  A2r,A2i;  /*  equation  (xxx)  in  the  text  */ 
float  B0,B1,B2;  /*  equation  (zsz)  in  the  text  */ 
float  wr [MAXSOX] , vi (MAXBOX) ;  /*  real  and  imag  upwash  at  boxes  */ 
float  Dr  {MAXDIM] ,  Di  (MAXDIM]  /  /*  The  real  and  imag  AIC  matrix  */ 
float  liftr,lifti;  /*  The  complex  lift  in  cartesian  form  */ 
float  liftm, liftp;  /*  The  complex  lift  in  polar  form  */ 
struct  element  box  [MAXJBOX] ;  /*  box  geometric  data  */ 
struct  trapezoid  wing;  /*  wing  geometric  data  */ 

seconds  a  PAUSE^ON_OUX?UT ;  /*»  seconds  the  program  will  pause  */ 

if  ( (odat^fopea ("dl, TRASH", "w") ) NULL) 

( 

print f  (*‘\ncannot  open  file  di .  TRASH  for  output \n")  / 
exit  (0) ; 

) 

printf (*\n\n  Auxiliary  data  placed  in  file  { dl. TRASH ) \n") ; 

if  ( (aicdat.«f  open  (Ml. AIC",  "w") )  mm  m) 

( 

printf ("\ncannot  open  file  dl.AlC  for  output \n" ) ; 
exit  (S‘,  ; 

) 

printf  (*\n  AIC  placed  in  file  (dl .AIC] \n") ; 
printf  P\n  HAXBOX:  %d  MAXDIM:  %d\n%  MAXBQX,  MAXDIM) ; 

/*  BEGIN  INPUT  */ 

printf  ("“SnBegin  input \n") ; 

ierr  «  read_input (odat , 4M, 4k, fcb, twing) ; 


no 
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px'intf  ("Input  complete\n") ; 

ierr  =  quadrilateral (odat, wing. xible, wing. yible, 

wing . xibte , wing . y ibte , 
wing . xobte , wing . y obte , 
wing . xoble , wing . y oble , 
fifactor, SxO, &y0) ; 

wing. area  =  factor/ 

printf ("\nWing  area  used  to  non-dimensionalize  lift  is:%12.4e\n", 
wing. area) ; 

printf ("The  wing  centroid  is  at  x:  %f  and  y:  %f\n",xO,yO) ; 
sleep (seconds) ; 

/*  INPUT  IS  NOW  COMPLETE  ,  DISCRETIZE  THE  WING  INTO  BOXES  */ 

printf ("Begin  discretizing  the  wing.\n*); 

ierr  =  discretize (odat, wing, box) ; 

printf ("Discretization  is  now  complete. \nw) / 

/*  non-dimensional  variables  are  computed  */ 

beta2  »  1-M*M/ 
b2  =  b*b; 

printf ("  reduced  freq  (k) :  %12.4e  \nw,k); 
printf ("  betaA2  :  %12.4e  \n*,beta2) ; 

/*  Zero  out  the  [D]  matrix  of  AIC  coefficients  */ 

for  (rb<aO ;  rb<wing .  total^boxes;  -M-rb) 

{ 

for ( sb»0 1 sb<wing . total_boxes / ++sb) 

( 

vector__index  ®  rb*wing.total_boxes+ab/ 

Dr [vector^index]  ®  0.0/ 

Di  (vector__index]  »  0.0/ 

) 

} 

ic  *>  0/ 
printf ("Xn") / 

I*  PROCEED  TO  FORMULATE  THE  COMPLEX  AIC  MATRIX  D!I,JJ  */ 

for  (symm»0 ;symm<®ABS (wing. syaa)  ;++symm)  /*  consider  symmetry  */ 

{ 

for  (rb=*G  ;  rb<wing .  total_boxes  /  ++rb)  /*  receiving  box  index  */ 

i 


ill 
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for  (sb=0;sb<wing.total__boxes;++sb)  /*  sending  box  index  */ 

{ 

/*  Kbar  is  defined  as  equation  285  */ 

/*  compute  Kbar  for  left  terminus  of  doublet  segment  */ 

xO  =  box[rb] .xc  -  box[sb].xi; 

if (symm==0) yO  =  box[rb].yc  -  box[sb] .yi; 

else  yO  =  box[rb] .yc  +  box[sb] .yo; 

fprintf (odat, 

"\nLEFT  symm:  %2d  rb:  %2d  sb:  %2d  xO:  %12.4e  yO:  %12.4e\n", 
symm,rb,  sb,xO,yO)  / 

ierr  =  Kbar (odat, M,k, xO,yO, &KbarrL, SKbariL) ; 

/*  compute  Kbar  for  midpoint  of  doublet  segment  */ 
xO  =  box[rb].xc  -  box[sb] .xm; 
if (symm==0) yO  =  box[rb] .yc  -  box[sb] .ym; 
else  yO  ~  box[rb].yc  +  box[sb].ym; 
fprintf (odat, 

"\nCENTER  symm:  %2d  rb:  %2d  sb:  %2d  xO:  %12.4e  y0:%12. 4e\n", 
symm, rb,  sb,xO,yO) ; 

ierr  ■  Kbar (odat, M,k, xO, yO, &KbarrC, &KbariC) ; 

/*  compute  Kbar  for  right  terminus  of  doublet  segment  */ 

xO  »  box[rb].xc  -  box[sb].xo/ 

if  (syrora^O)  yO  **  box[rb].yc  -  box[sb].yo; 

else  yO  =  box[rb] .yc  +  box[sb].yi; 

fprintf (odat, 

"\nRIGHT  symm:  %2d  rb:  %2d  sb:  %2d  xO:  %12.4e  y0:%12.4e\n", 
symm, rb, sb,xO,yO) ; 

ierr  =*  Kbar (odat, M,k, xO,yO, &KbarrR, SKbariR) ; 

dLx  «  box[sb] .xc-box[sb] .xi;  dLy  ■  box [sb] .yo-box[sb] .yi/ 

L  =  sqrt ( (double) (dLx*dLx+dLy*dLy) ) /2.0; 

L2  ■  L*L; 

/*  set  the  sweep  angle  here  for  the  doublet  line  segment  */ 
if(aymm®=0)  si  =  (box [sb] . yo-box [sb] ,ym) /L; 
else  si  “  (box[sb] .ym-box[sb] .yi) /L;  /*  left  half  of  wing  */ 
sl2  ■  sl*sl; 

*13  *  sl*&12/ 

if(symm««0)  Yt  »  box[rb].yc  -  box[ab].ym/  /*  local  y  axis  */ 
else  Yt  =  box[rb].yc  +  box[sb].ym/  /*  left  half  of  the  wing  */ 
Yt2  =  Yt*Yt; 

/*  The  real  and  imaginary  components  of  the  A  coeffiecients 
of  equations  288-290:  */ 

AOr  ■  KbarrC; 

AOi  *  KbariC/ 
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Air  =  (KbarrR-KbarrL) / (2*L) / 

Ali  =  (KbariR-KbariL) / (2*L) ; 

A2r  =  (  KbarrL  -  2.0*KbarrC  +  KbarrR  )/(2.C*L*L); 

A2i  =  {  KbariL  -  2.0*KbariC  +  KbariR  )/{2.0*L*L); 

factO  =  (sl2*L2-2.0*Yt*sl*L+Yt2) / (sl2*L2+2 . 0*Yt*sl*L+Yt2) / 
factO  =  log ((double) (factO));  /*  nat  log  */ 

BO  =  (2 . 0*L) / ( Yt2-I2*sl2) ;  /*  equation  296)  */ 
factl  =  (0.5/sl2) *factO; 
fact2  =  (Yt/ si) *B0  j 

B1  =  (factl+factl) /  /*  equation  298  */ 
factl  =  2.0*L/sl2; 
fact2  =  (Yt/sl3)  *£actf»; 
fact 3  =  (Yt2/sl2) *B0; 

B2  =  (factl ffact2+fact3) /  /*  equation  300  */ 
factor  =  (-box [ab] . chord/ (8 . 0*PI) ) ; 

/* 

Here,  we  will  co  p  'ss  square  m&ir.ricies  Dr  and  Di  into  vectors. 
The  matricies  Dr  and  Di  are  the  real  and  imaginary  parts  of  [D] . 
<w]  =  [D](Cp] 

{w}  is  a  vector  of  non-dim  upwash  w/U  at  the  control  points. 

{Cp}  is  the  vector  of  non-dimensional  pressure  dP/ (rho*UA2) . 

[DJ  is  compressed  into  (Dr)  and  (Di)  in  row  packets. 

*/ 

vector_index  =  rb*wing. total_boxes+sb; 
if  (syram==30)  /*  for  any  symmetric  case  */ 

{ 

Dr [vector_index]  ®  Dr [vector_index]  + 

factor*  (B0*AQr4-Bl*Alr'fB2*A2r)  / 

Di  (vector__indexJ  «  Di  [vector__index]  + 

factor*  (£G-A0i-*-5i '‘Aii+Bz  wA2i )  ; 

} 

else  if  (symm«*“l& &wing. symm>0)  /*  left  half  wing  is  aymm  */ 

( 

Dr  [vector_index]  **  Dr [vector_index]  + 

factor* (B0*A0r+Bl*Alr+B2*A2r) / 

Di  [vector__index]  «  Di  [vector_index]  -f 

factor* (B0*A0i+Bl*Ali+B2*A2i) ; 

) 

else  if  (symm®*=l& Swing. symm<0)  /*  left  wing  is  anti-symm  */ 

< 

Dr [vector_index]  ■  Dr [vector_index]  - 

factor* (B0*A0r+Bl*Alr+B2*A2r) ; 

Di [vector_index]  »  Di [vector_indox]  - 

factor* (B0*A0i+Bl*Ali+B2*A2i) ; 
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} 

else 

{ 

printf ("\n  CONFUSED  about  the  wing  symmetry...  \nexit\n") / 
exit (0) ; 

} 

}  /*  end  of  loop  on  sb  */ 

if (++ic>10) 

{ 

printf ("\n") ; 

'.0=0; 

} 

printf  ("  %d"’,rb); 

}  /*  end  of  loop  on  rb  */ 

}  /*  end  of  loop  on  symm  */ 
printf  ("\n\n'r) ; 

/*  Print  out  the  [D]  matrix  of  AIC  coefficients  */ 
for ( rb=0 ; rb<wing . total  Jboxes ; ++rb) 

{ 

for (sb«0 / sb<wing . total  Jboxes ; ++ab) 

{ 

vector_index  »  rb*wing. total  J>oxes+sb; 

fprintf  (aicdat, "row:  %5d  col:  %5d  ind:  %10d  %15.7e  %15.7e\n*rf 
rfo, sb, vector_index,  Dr [vector_index] , Di [v«ctor_ihdex] ) ; 

} 

> 

/*  Input  the  {w>  boundary  condition,  return  #  of  monomials  */ 
if (  (ierr  «  bc(odat,k,wr,wi,box, wing) )  *»  0  ) 

{ 

printf ("\nNo  upwash  specified  and  no  pressure  computed\n*) ; 
f close (odat) ; 
f close (aicdat) ; 
exit (0) ; 

) 

else 

{ 

printf  (*\n  Upwash  specified  and  pressure  will  be  computedXn") ; 

) 

for  (i»0;i<ving.total_boxes;'M-i) 

{ 

printf ("%5d  Real (w] :  %12.4e  Imaglv’:  %12.4e  \n" , 
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i,wr[i]  ,wi  [i] ) ; 

fprintf (odat, "%5d  Real[w]:  %12.4e  Imag[w] :  %l2.4e  \n ", 
i,wr [i]  ,wi  [i] ) ; 

} 

/*  solve  the  complex  problem  {w}  =  [D] {p}  */ 
printf ("\nSolve  the  complex  problem  {w}  «  [D]{p}\n")/ 
ierr=complex_solve (Dr,Di,wr, wi,wing.total_boxes) ; 
if  (lerr  1=0) 

{ 

printf ("\nerr or  number  %d  in  complex_solve\n", ierr) ; 
exit ( 0 ) ; 

} 

/*  note:  the  [D]  matrix  is  now  the  complex  identity  matrix  */ 

/*  Print  out  the  pressrre  coefficients  and  sum  the  lift*/ 

printf ("\n\nCp  PRESSURE  COEFFICIENTS  (P=0 . 5*rho*UA2*Cp) n) / 

printf  (  "\nbox  ♦  (rea1  Cp)  (imag  Cp)  (box  area) \nw) ; 

f printf  (odat ,  w\n\.iCp  PRESSURE  COEFFICIENTS  (P«0.5*rho*UA2*Cp)  *)  ; 

f printf (odau, *\nbox  #  (real  Cp)  (imrg  Cp)  (box  area) Nn*) / 

liftr  *=  0.0/  lifti  =  0.0/ 

for  ( i- =0  /  i<wing .  total_boxeo  /  ++i ) 

{ 

printf (  "%5d  %12.4e  %12.4e  %12.4e\n*, 
i,wr [ij ,wi [i] ,b2*box[i] .area) / 
fprintf (odat,  w%5d  %12.4e  %3'>.^e  %12.4e\n*, 
i, wr  [i]  ,wi  (i]  ,b2*box [i)  .ar-a)  ♦ 
liftr  wr [i] *box [i] .area/ 

lifti  +«  wi[i.] *box(i]  .area/ 

> 

liftr  /*  wing.&x'ea / 
lifti  lv  wing. area/ 

ierr  polar (liftr, lifti, filiftm, iliftp) / 
liftp  *•  (180 . 0/PI } / 

printf ("\nThe  characteristic  length  b:  ",b)/ 
printf ("The  wing  area  used  to  non -dim  lift  is:%12.4e\n", 
wing. area) / 

printf ("Lift  Coefficient  *  C_L*q*A\n") / 

printf ("THE  COMPLEX  WING  LIFT  COEFFICIENT  (C_L)  IS:  "); 

printf ("( (%11.4e)  +  (%11. 4e)i)\n*, liftr, lifti) / 

printf  ("MAGNITUDE:  (%11.4e)  PHASE;  (%S.4f) deg\n'r,liftm,liftp)  / 


fclose (odat) / 
f close (aicdat) / 
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exit  (0)  / 
} 


FUNCTION :  r«ad_input ( ) 

♦include  "dl. define" 

♦include  <math.h> 

♦include  <stdio.h> 

♦include  <string.h> 

♦include  "dl . structure" 

r ead_input ( odat , M, k , b , wing ) 
float  *M,  *k,  *b; 
struct  trapezoid  *wing; 

FILE  *odat; 

{ 

FILE  *fopen(); 

FILE  *idat/ 
char  csymm; 
char  line [200]; 
int  ierr; 
int  n; 
float  X/  y; 

if ( (idat-fopen ("dl . INPUT", "r") ) —  NULL) 

{ 

print f ("\ncannot  open  file  dl. INPUT  for  input\n") ; 
exit (0) ; 

} 

printf ("Input  data  will  be  read  from  file  [dl. INPUT] \n") ; 

/*  BEGIN  INPUT  */ 

if(  in_line  (idat,  line)  “0  )  /*  read  title  line  */ 

{ 

fprintf (odat, "\ntext ;  [%a] ",  line) ; 
printf ("TITLE :\ntextj  [%a] \n", line) / 

} 

if (  in-line (idat, line) *»0  )  /*  read  characteristic  length  */ 
( 

sscanf (line, "%f", tx) ; 

*b  ■  x; 

printf ("characteristic  length:  %f\n",*b); 
fprintf (odat, "characteristic  length:  %f\a",*b); 

) 

if (  in  line (idat, line) )  /*  read  Mach  number  */ 
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{ 

sscanf (line, "%f",  &x)  ; 

*M  =  x; 

printf ("Mach:  %-p\n",*M); 
fprintf  (odat,  "Mach:  %f\n",  *M) ; 

} 

if  (  in__line  (idat,  line)  =0  )  /*  read  reduced  frequency  */ 

{ 

sscanf (line, "%f", &x) ; 

*k  =  x; 

printf ("reduced  frequency:  %f\n",*k); 
fprintf (odat, "reduced  frequency:  %f\n",*k); 

} 

if  (  in-line (idat, line) *=0  )  /*  read  symmetry  flag  -1,  0  or  +1  */ 
{ 

sscanf (line, "%c", ficsymm) ; 
line[0]  =  csymm;  line[l]  =  '\0'; 
wing->symm  =0; 
if  (  strcmpdine,  "s")“»0  ) 

{ 

wing~>symm  =  1/ 

printf ("Assume  symmetry  about  the  x  axis\n") / 
fprintf (odat, "Assume  symmetry  about  the  x  axis\n") ; 

} 

else  if (  strcmp(line,"a")**0  ) 

( 

wing->syitun  «  -1; 

printf ("Assume  anti -symmetry  about  the  x  axis\n") ; 
fprintf (odat, "Assume  anti-symmetry  about  the  x  axie\n")/ 

} 

else 

{ 

wing->symm  *  0; 

printf ("Assume  no  symmetry  about  the  x  axis\n") ; 
fprintf  (odat,  "Assume  no  symmetry  about  the  x  axis\n")  / 

> 

} 

if(  in  line (idat, line) «*0  )  /*  read  inbrd  lead  edge  coord  */ 

{ 

sscanf (line, "%f  %f",4x,4y); 

%ring->xible  *  x; 
wing->yible  ■  y; 

printf ("inboard  leading  edge:  x  %f  y  %f\n", 
wing->xible, wing->yible) ? 
fprintf (odat, "inboard  leading  edge:  x  %f  y  %f\n", 
wing->xible, wing->yible) / 
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wing->xible  /=  *b; 
wing->yible  /*=  *b; 

} 

if  (  in_line  (idat ,  line)  =»«0  )  /*  read  inbrd  trail  edge  coord  */ 

{ 

sscanf (line, "%f  %f",  &X,  &y)  ; 
wing->xibte  =  x; 
wing->yibte  =  y; 

printf ("inboard  trailing  edge:  x  %f  y  %f\n", 
wing->xibte, wing->yibte) / 

fprintf (odat, " inboard  trailing  edge:  x  %f  y  %f\n", 
wing->xibte, wing->yibte) ; 
wing->xibte  /=  *b; 
wing->yibte  /=  *b; 

} 

if  (  in_line (idat, line) =®0  )  /*  read  outbrd  trail  edge  coord  */ 

{ 

sscanf (line, "%f  %f",&x,&y); 
wing->xobte  =  x; 
wing->yobte  «=  y; 

print f ("outboard  trailing  edge:  x  %f  y  %f\n", 
wing- >xobte, wing- >yobte) ; 

fprintf (odat, "outboard  trailing  edge:  x  %f  y  %f\n", 
wing- >xobte, wing- >yobte) ; 
wing->xobte  /«=  *b; 
wing->yobte  /**  *b; 

) 

if  (  in_line  (idat,  line)  »»0  )  /*  read  outbrd  lead  edge  coord  */ 

{ 

sscanf (line, "%f  %f",4x,6y); 
wing->xoble  ■  x; 
wing->yoble  «  y; 

printf ("outboard  leading  edge:  x  %f  y  %f\n", 
wing->xoble,wing->yoble) ; 

fprintf (odat, "outboard  leading  edge:  x  %f  y  %f\n*, 
wing->xoble, wing->yoble) ; 
wing->xoble  /*  *b; 
wing->yoble  /•  *b; 

) 

if  (  in-line (idat, line) -»0  )  /*  read  number  of  boxes  in  x  */ 

t 

sscanf (line, "%d",  4n) / 
wing->num_box_x  «=  n? 

printf ( "number  of  boxes  in  the  x  direction:  %d\n", 
ving->nua_box__x) ; 

fprintf (odat, "number  of  boxes  in  the  x  direction:  %d\n", 
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wing-  >num_box__x ) ; 
if  (wing->num_box_x  >  MAX_DIV_X) 

{ 

print f (v\n  EXCEEDED  MAXIMUM  DIMENSION  ON  DIVISIONS  IN  X\nw) ; 
exit  (0) / 

} 

} 

if  (  in_line  (idat,  line)  =0  )  /*  read  number  of  boxes  in  y  */ 

{ 

sscanf (line, "%d", &n) ; 
wing- >num_box_y  =  n; 

print f ("number  of  boxes  in  the  y  direction:  %d\n", 
wing->num_box_y) / 

fprintf (odat, "number  of  boxes  in  the  y  direction:  %d\n", 
wing-  >num_box_y )  / 
if (wing->num_box_y  >  MAX  DIV_Y) 

{  ’ 

print f  ("\n  EXCEEDED  MAXIMUM  DIMENSION  ON  DIVISIONS  IN  Y\n") ; 
exit  ( 0 ) ; 

} 

} 

wing->totsi_boxea  =  ving->num_box_x  *  wing- >num_box_y ; 
if (wing->totalJboxes>MAXBOX) 

( 

print f ("\n  EXCEEDED  MAXIMUM  DIMENSION  ON  BOXES \n") ; 
exit  (0) ; 

> 

f close (idat) / 

/*  INPUT  IS  NOW  COMPLETE  */ 

return ( 0 ) ; 

) 


EBMfiTIOH;. 


/**************************************************************** 
‘Calculate  the  coordinates  of  the  discretised  wing.  * 

********************************** »*****************************/ 

♦include  "dl. define" 

♦include  <math.h> 

♦include  <stdio.h> 

♦include  "dl . structure" 
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int  di s cr ©t i z  e ( odat , wing , box ) 

FILE  *odat; 

struct  element  box [ ] ; 

struct  trapezoid  wing; 

{ 

int  quadrilateral () / 
int  ierr; 

int  rx,ry;  /*  x  index  and  y  index  */ 
int  box_index; 

float  dxib, dxob;  /*  delta  x  on  the  inboard  and  outboard  chords  */ 

float  dy;  /*  delta  y  is  const ant  over  the  span  */ 

float  xl,x2,x3,x4; 

float  yl,y2,y3,y4; 

float  ix, iy; 

float  area, total_area; 

total_area  “  0,0; 

/*  delta  x  along  the  inboard  and  outboard:  */ 

dxib  *  (wing. xibte-wing. xible) / (float) wing.num_box_x; 

dxob  *  (wing.xobte-wing.xoble)  /  (float) wing. num_box_x; 

/*  delta  y  of  all  boxes:  */ 

dy  ®  (wing. yob le-ving.yible) / (float) wing, numjbox^y; 
f print f  (odat, ''dxib:  %f  dxob:  %f  dy:  %f\n",dxib,dxob,dy) ; 

f or ( ry«Q ; ryewing . num  box_y;++ry)  /*  box  y  index  */ 

{ 

for  (rx»0;rx<wing.namjbox_jK;+*rx)  /*  box  x  index  */ 

{ 

box^index  *=  ry *wing .  num^box^x  +  rx; 

/*  printf ("rx:  %3d  ry;  %3d  receive^index:  %4d\n", 
rx, ry , box_index) ;  *  / 

/*  inboard  leading  edge  coordinates  of  receiving  box:  */ 
ix  *  (float) rx; 

iy  *  (float)  ry/ (float)  wing. num_box__y; 
xl  «  ix*dxib  +  ix* (dxob-dxib) *iy;  /*ible*/ 
yl  -  ry*dy; 

/*  inboard  trailing  edge  coordinates  of  receiving  box:  */ 

ix  «  (float)  (rx-fl) ; 

iy  *•  (float) ry/ (float) wing. num_box_y; 

x2  «  ix*dxib  +  ix* (dxob- dxib) *iy?  /*ible*/ 

y2  *  ry*dy; 

/*  outboard  trailing  edge  coordinates  of  receiving  box:  */ 
ix  **  (float)  (rx+1); 

iy  *  (float)  (ry+1)  /  (float)  wing.  numj5ox__y; 
x3  •  ix*dxib  +  ix* (dxob- dxib; *iy;  /*ible*/ 
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y3  =  (ry+1)  *dy; 

/*  outboard  leading  edge  coordinates  of  receiving  box:  */ 
ix  =  (float) rx; 

iy  =  (float)  (ry+1)  /  (float)  wing. num__box_y; 
x4  «  ix*dxib  +  ix* (dxob-dxib) *iy;  /*ible*/ 
y4  =  (ry+l)*dy; 

/*  coord  of  the  receiving  control  pt  at  3/4  chord  centerspan  */ 
box [box_index] .xc  =  (  xl+0 , 75* (x2-xl)  +  X4+0.75* (x3-x4)  )/2.0; 
box [box_index] . yc  =  (y2+y3) /2 . 0 ; 

/*  inboard  coord  of  the  sending  doublet  line  along  1/4  chord  */ 
box [box_index] .xi  =  xl+0.25* (x2-xl) ; 
box[box_index] .yi  =  yl; 

/*  outboard  coord  of  the  sending  doublet  line  along  1/4  chord  */ 
box[box_index] .xo  =  x4+0.25*  (x3-x4)  ; 
box (box_index ] .yo  »  y4; 

/*  midspan  coord  of  the  sending  doublet  line  along  1/4  chord  */ 
box [box_index] . xm  =  (box [box_index] . xi+box [box__index] . xo) /2 . 0 ; 
box [box_index] . ym  =  (box [box_index] . yi+box [box_index] . yo) /2 . 0; 
/*  average  chord  */ 

box [box_index]  . chord  =  ( (x2-xl) + (x3-x4) ) /2. 0; 

/*  box  area  and  x  and  y  coordinates  of  the  box  centroid  */ 
ierr  =  quadrilateral (odat,xl,yl,x2,y2,x3,y3,x4,y4, 

&araa, Six, &iy) ; 
box [box_index] .area  «  area; 
box [box_index] .xcent  “  ix; 
box[box_index] .ycent  “  iy; 
total_area  +■  area; 

fprintf (odat, w\n  BOX:  fcSdXn", box__index) ; 

fprintf (odat, "  xl:  %8.4f  x2:  %8.4f  x3;  %8.4f  x4:  %8.4f\n*, 
xl , x2 , x3 , x4 ) ; 

fprintf (odat, "  yl:  %8.4f  y2:  %8.4f  y3:  %8.4f  y4:  %8.4f\n", 
yl,y2,y3,y4) ; 

fprintf (odat, "  3/4  chord  midspan  x:  %f  y?  %fXn*% 
boxtbox_index] . xc,box[box_index) .yc) ; 
fprintf (odat, *  1/4  chord  inboard  x:  %f  y:  %f  \n", 
box  (box__index)  .xi,box(box_jLndex]  .yi) ; 
fprintf (odat, "  1/4  chord  midspan  x:  %f  y:  %f  Xn*, 
box [box_index) . xm, box [box_index] . ym) ; 
fprintf  (odat, "  1/4  chord  outboard  x:  %f  y:  %f  \n"’, 
box [box_index] . xo,box [box_index] .yo) * 
fprintf  (odat, "  average  chord:  %f  Xu*, box [box__xndex]  .chord) ; 
fprintf (odat, *  area:  %f  \n*,area); 

fprintf  (odat, 0  x  centroid:  %f  y  centroid:  %f\n*,ix,iy) ; 

)  /*  end  of  loop  on  rx  */ 

}  /*  end  of  loop  on  ry  */ 
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fprintf  (odat,  "\nTOTAL  AREA  OF  WING:  %f\n*,total_area) ; 
return { 0 ) ; 

} 


FUNCTION :  KbarO 

/********************************************************!>******* 

*  This  subroutine  computes  K  bar.  * 

*  K  bar  is  given  as  equation  273  in  the  text) .  * 

*  K1  is  given  as  equation  266  in  the  text.  * 

**************************************************************** I 

♦include  Ml. define" 

♦include  <math.h> 

♦include  <stdio.h> 

int  Kbar (odat,M,k,xO,yO, Kbr,Kbi) 

FILE  *odat; 

float  M; /*  mach  *J 

float  k;/*  reduced  freq  wb/U  */ 

float  xO,yO;/*  non-dimeneional  (x-xi) /b  and  (y-eta) /b  */ 
float  *Kbr,*Kbi;/*  return  these  values  */ 

{ 

int  ierr,Il(); 
float  ul,kl; 
float  alpha, bet a2; 
float  Kir, Kli, factor; 
float  exr,exi,eur,eui; 
float  Hr,  Hi; 
float  Ul; 

fprintf  (odat,  "In  Kbar  no%r\n")  / 
fprintf  (odat,  *M:  %12 , 4e\n*',M) ; 
fprintf (odat, *k:  %12 . 4e\n*, k) ; 
fprintf  (odat,  "xO:  %12.4e  yO:  %12.4e\n",x0,y0) ; 

beta2  = 

kl«=  k*ABS  (yO) ; 

fprintf (odat, *kl:  %12.4e\n*,kl) ; 

if  ( (ABS (yO) )<EPS)  /*  if  yO  ■  aero,  ve  nesd  to  take  care  for  ul  */ 

( 

if (x0>0. 0) (Ul-EIGM; } 
else  (ul-BIGP;) 

) 


FUNCTION:  110 


else  /*  yO  is  not  equal  to  zero  */ 

{ 

ul  =  (M*sqrt ( (double)  (x0*x0+beta2*y0*y0) ) -xO) / (ABS (yO) *beta2) ; 

} 

alpha=  ~xO*k; 

exr=  cos { (double) (alpha) ) ; 
exi=  sin ( (double) (alpha) ) ; 

ierr  =  II (odat,ul,kl, &IIr, &Ili) ;  /*  compute  the  II  integral  */ 

/*  compute  K1  =  (Kir  +  i*Kli)  */ 
if (  ul>=BIGP  | |  ul<=BIGM  ) 

{ 

fprintf  (odat,  "’BIG  ul:  %12. 4e\n", ul) ; 

Kir  =  -Hr; 

Kli  =  -Ili; 

} 

else  if (  ul<BIGP  &&  ul>BIGM  ) 

{ 

fprintf  (odat,  ''bounded  ul:  %12.4e\n",ul) ; 

alpha  =>  -kl*ul; 

eur  ®  cos( (double) (alpha) ) ; 

eui  s  sin ((double) (alpha)); 

factor  =  (-1.0) * 

(  M*ABS (yO) /sqrt ( (double) ( (x0*x0) «• (beta2*y0*y0) ) )  ) 

/ 3 qrt ( (double) (1.0+ul*ul) ) ; 

Kir  »  factor*eur  -  Hr; 

Kli  »  factor*eui  -  Ili; 

> 

else 

( 

printf (*\aConfused  with  ul  in  function  KbarVa*)  / 

} 

/*  compute  Kbar  »  (Kbr  +  i*Kbi)  */ 
ierr  *  cmult (Klr,Kli,exr,exi,Kbr,Kbi) ; 
return (ierr) ; 

) 


FUNCTION :  IK) 

/A******************************************************* 

*  This  soubroutine  computes  the  integral  II.  * 

*  II  is  given  as  equation  267  in  the  text  * 

*  II  computes  the  II  integral  for  any  ul.  * 

************* ******************* ****************»*******/ 
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FUNCTION:  I1Q 


♦include  <math.h> 

♦incltide  <stdio.h> 

♦include  Ml.  define* 

int  Il(odat,ul,kl,Ilr,Ili) 

FILE  *odat; 

float  ul,kl;/*  input  these  values  */ 
float  *Ilr, *Ili;/*  return  these  values  */ 

{ 

int  ierr,  cmult(); 
float  Ul,  I real; 

fprintf  (odat,  *  (IN  II)  ul:  %12.4e  kit  %12.4o\n*,ul,kl) / 
if (ul>®0.0) 

{ 

if (ul>«BIGP)  /*  +BIG  <  ul  */ 

{ 

*Ilr  -  0.0; 

*IXi  *  0.0; 

) 

else  /*  0  <•  ul  <  +BIG  */ 

{ 

ierr  -  il (odat,ul,kl,Ilr, Ili) / 

} 

) 

else  /*  ul  <  0  */ 

< 

if(ul<«BIGM)  /*  ul  <  -BIG  */ 

{ 

fprintf (odat,*ul  <  -BIG  \»*) ; 

Ul  -  0.0; 

ierr  -  il  (odat, 01, kl.  Hr,  Hi) ; 

*Ilr  *«  2.0; 

*Ili  2.0; 

) 

else  /*  -BIG  <  ul  <  0  */ 

( 

/*  coapute  II  for  0<u<i»finity:  */ 
fprintf (odat,*  -BIG  <  ul  <  0\n*) ; 

Ul  -  0.0; 

ierr  «•  il  (odat,Ul,kl,  Hr,  Hi) ; 

Ireal  -  *Ilr; 

Ul  «  (-1.0) *ul; 

ierr  »  il  (odat,Ul,kl,  Ilr,  Ili) ; 
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FUNCTION  il() 


*Ilr  =  (2. 0*lreal)  -  (*Ilr) / 

} 

} 

return (ierr)  ; 

} 


FUNCTION  il Q 

/****************************$*********************************** 
♦Function  to  compute  the  II  integral  for  ul  >»  0* 

a***************************************************************/ 

int  il (odat,ul,kl,ilr,ili) 

FILE  *odat; 
float  ul,kl; 
float  *ilr,*ili; 

{ 

int  ierr,  cmult  () ,  Jl  () ; 
float  ir, ii, alpha, factor, er,ei; 
float  jlr, jli/ 

fprintf (odat, " (IN  il)  ul:  %l2.4e  kl:  %12.4e\n",ul,kl) ; 
if  (uKO.O) 

{ 

printf(w\n  ul  cannot  be  less  than  zero...ul  »  %12.4e\n*,ul) ; 
exit (0) ; 

} 

ierr  **  Jl (odat,ul, kl, & jlr, 4 jli) ; 

factor  »  1.0  -  ul/sqrt  ( (double)  (1.0+ul*ui) ) ; 

ir  *  factor  +  kl*jli; 

ii  «=>  -kl*jlr; 

alpha  «  -kl*ul; 

er  *  cos ( (double) (alpha)); 

ei  *  ain( (double) (alpha)); 

ierr  -  cmult (er,ei,ir,ii,ilr, ili) ; 

return (ierr) ; 

) 


FUNCTION  JIM 

/*«***  ft  ft******************  «;«********»  **********  ********** 

*  This  subroutine  computes  the  integral  Jl.  * 
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FUNCTION  J10 


*  The  integrand  of  J1  is  appro-.  mated  as  an  * 

*  algebraic  expression  and  the*,  integrated.  * 

*  The  series  expression  for  J1  is  given  as  * 

*  equation  xxx  in  the  text.  * 

********************************************************  [ 

♦include  <math.h> 

♦include  <stdio.h> 

int  J1 (odat,ul,kl,  Jlr, Jli) 

FILE  *odat, 

float  ul,kl;/*  input  thesa  values  */ 
float  *Jlr,*Jli;/*  return  these  values  */ 

{ 

int  i,ierr; 

static  int  n  *  11; 

fxoat  j,zero, sjr, sji, jr, ji; 

double  djr, dji; 

static  float  c  =  0.3'72; 

static  float  a[ll]  *  {  0.24186198, 

-2.7918027, 

24.991079, 

-111.59196, 

271,43549, 

-305.75288, 

-41.183630, 

545.98537, 

-644.78155, 

328.72755, 

-64.279511  }/ 

fprintf  (odat, "  (IN  Jl)  ul:  %12.4e  kl:  %12. 4e\n*,ul,kl) ; 

ierr  «  0; 

zero  *=  0.0; 

djr  *  (double) zero; 

dji  «  (double) zero; 

for ( i*l ; i<»n ; ++i ) 

( 

j  «  a [i-l]*exp( (double) (-i*c*ul) ) / 

((float) (i*i)*(c*c)+(kl*kl))/ 
sjr  *  j*  (fl oat) i*c; 
sji  =  -j*kl; 

djr  »  djr  +  (double) sjr; 
dji  ®  dji  +  (double) sji; 
fprintf (odat, "a (%2d) ;  %12.4e\n", i, a[i-lj ) ; 
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FUNCTION  bc() 


} 

*Jlr  :*  ( float )djr/ 
*Jli  -  (float) dji; 

return (ierr) ; 

} 


FUNCTION  beO 

♦include  MI.  define" 

♦include  <math.h> 

♦include  <stdio.h> 

♦include  ^string. h> 

♦include  "dl . structure" 

int  he (odat,k, wr,wi, box, wing) 

FILE  *odat; 
float  k; 

float  wr [ ] ,  wi  t ] ; 
struct  element  box [ ] ; 
struct  trapezoid  wing; 

{ 

FILE  *fopen(); 

FILE  *idat; 

char  c, flag [2 j , line [200]; 
int  ierr, do_flag, number; 
int  i,j,px,py; 
float  a,  x,  y,  sumr,  sumi; 
float  power (); 

struct  polynomial  poly  (MAXJ?0L¥) ; 

if  ( (idat-fooen  ("be.  INPUT",  "r") ) «»  NULL) 

{ 

print f ("\noannot  open  file  be. INPUT  for  input \n") ; 
exit  (0) ; 

> 

printf ("Boundary  data  will  be  read  from  file  [dl. INPUT) \n") ; 
printf ("reduced  frequency  (k) :  %12.4e\n",k) ; 

number  «  0; 

/*  BEGIN  INPUT  */ 

do{ 

do_flag  »  1; 
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FUNCTION  bcQ 


if  (  in__line  (idat, line)  «»0  ) 

{ 

/*  fprintf (odat, "\ntext:  [%s] ",  line) /  *J 
l*  printf ("text:  [%s] \n", line) ;  */ 
sscanf (line, "%e", sc) ; 
flag  [0]  **c; 
flag[l]«'\0'; 

if(  strcmp(flag, } /*  ©nd  of  data  */ 

{ 

do_flag*0 ; 

} 

if(  strcmp(flag,"’l#')«»0  ) /*  a  line  of  data  is  in  line[]  */ 

{ 

sscanf  (line,  "%c  %f  %d  %d",6c, &a,Spx,£py)  / 
poly [number]  .a  *  a;  /*  coefficient  */ 
poly [number] .px  »  px;  /*  power  of  x  *J 
poly [number] .py  -  py;  /*  power  of  y  */ 
printf("  [%c]  a:  %f  px;  %d  py:  %d\n*, c,a,px,py) ; 

++number; 
do  flag  »  1/ 

) 

} 

)while(do_>flag»»l) ; 
f close (idat) ; 

if (number««0) return  (number) ;  /*  return  sero  if  no  data  input  * / 

/*  Compute  the  upwaah  at  each  box  3/4  chord  */ 

print f("\n  compute  the  upwash  at  %d  control  points' \n", 
wing.totaljboxes) ; 
for (i*0;i<wing. total  boxes ;++i) 

< 

x  »  box[i] .xc; 
y  -  box[i] .yc; 
sumr  “  0.0; 
sum!  •  0.0; 
f or < i«0 ; j<number ; +♦ j) 

{ 

sumr  **  sumr  -f  (poly [j]. a  *  rolytjJ-P*)  * 

power (x, (poly [j] .px-l))  *  power (y, poly [j] .py) ; 
sum!  *  aumi  +  k  *  polytj] .a  * 

power (x, poly [j] .px)  *  power (y, poly (j] .py) ; 

> 

wr[i]  «  sumr; 
wi[i]  ■  sumi; 
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FUNCTION  complexjsolveO 


fprintf (odat, "  %5d  Real[w]:  %l2.4e  Imag[w] :  %12.4e\n"/ 
i,wr[i] / wi [i] ) ; 

} 

return (number) ; 

} 


FUNCTION  complex^aolve ( ) 

/* 

solve  is  a  function  to  solve  the  complex  linear  system  [a] {x}={c> 
using  Gaussian  elimination  and  back  substitution  with  pivoting 
on  each  step,  [a]  is  input  in  vector  form  a(k)  *  a(i,j)  where 
k=i*nc+j.  nc  is  the  utilized  portion  of  a[nc] [nc]  and  c[nc) 
where  nc  <=  the  declared  dimension.  The  complex  solution  {x} 
is  returned  in  {c}. 

The  real  and  imaginary  parts  of  a[] []  are  designated  as 
ar[]  and  ai[].  Likewise  for  the  c  vector. 

*/ 

♦include  <stdio.h> 

♦define  ABS (x)  (((x)<0)  ?  -(x)  :  (x) ) 

int  complex_3olve (ar, ai, cr, ci, nc) 
int  nc; 

float  ar  [] ,  ai  [J  ,cr[]  ,ci[]  / 

{ 

int  ierr, i, j, i2,nr, ir, jc; 
int  cdiv ( ) ,  cmult ( ) ; 

float  rmaxabs ,  tmpabs,  tmpr ,  tmpi,  dr ,  di ,  big*l .  0e+20 ,  eps**l .  0e-20 / 
float  cmag ( ) ; 

/*  Consider  nc  rows:  */ 
for  (i«0 ;  i<nc;  ++i) 

{ 


/*  find  max  diagonal  value  and  switch  rows  */ 
nr«*i; 

rmaxabs»cmag(ar [i*nc+i] ,ai[i*nc+i] )  ; 
for (i2«i;i2<nc/++i2) 

{ 

tmpabs«rmaxabe-cmag(ar[i2*nc+i] ,ai [i2*nc+i] ) ; 
if  (tmpabs<0.0) 

{ 
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FUNCTION  complex_solve() 


nr*i2; 

rmaxabs^cmag  ( ar  f  i2  *nc+i  ] , ai [ i2  *nc+i  J  >  ; 

} 

}  /*  end  i2  loop*/ 

for ( j»0; j<nc;++ j) 

{ 

tmpr*ar [nr*nc+ j ] ; 
tmpi»*ai  [nr*nc+j] ; 
ar [nr*nc+j]*ar [i*nc+j] ; 
ai [nr*nc+j]»ai [i*nc+j] ; 
ar [ i *nc+ j ] *tmpr ; 
ai  [i*nc+j]«»tmpi; 

} 

tmpr=*cr  [nr] ; 
tmpi=ci  [nr] ; 
cr[nr]*cr [i] ; 
ci[nr]»ci [i] ; 
cr  [i]**tmpr; 
ci[i]*tmpi; 

/*  rows  have  been  switohed  */ 

dr-ar [i*nc+i] ; 
di»ai [i*nc+i] ; 

if  (ABS (cmag (dr, di) )<«epa) return (1) ; 
for ( j*i; j<nc;++ j) 

{ 

ierr«cdiv(ar [i*nc+j] , ai[i*nc+j] , dr,dif 4tapr,ttapi) / 
ar[i*nc+j]  «  tmpr; 
ai[i*nc+j]  «*  tmpi; 

) 

ierr»cdiv(cr  [i] ,ci  [i]  ,dr,di,  4tmpr,  fitmpi) ; 
cr[i]  «  tmpr; 
ci[i]  «  tmpi; 

if(  ABS (cmag (cr[i],ci[i)))  >•  big  ) return (5); 

for (ir»0;ir<nc;++ir) 

{ 

if (ir-*i) continue; 
dr  ■  ar[ir*nc4i]; 
di  -  ai[ir*nc+i]; 

if(  ABS (cmag (dr, di) )  >-  big) return (2) ; 
if(  ABS (cmag (cr [i] ,ci[i] ) )  >*  big) return <4 ) ; 
for ( jc*i ; jc<no; 4+ jc) 

< 

ierr-cmult  (dr,di,ar [i*nc4jc] ,  ai  [i*nc4 jc] ,  6 tmpr,  4 tmpi) 


FUNCTION  cacipiex^mathO 


ar [ir*nc+jc]  -=  tmpr; 
ai[ir*nc+jc]  -=  tmpi; 

}  /*  end  jc  loop  */ 

ierr=cmult (dr,di, cr[i] ,ci [i] , fitmpr, &tmpi) ; 
cr[ir]  -=  tmpr; 
ci[ir]  -=  tmpi; 

}  /*ir  loop*/ 

}  /*  i  loop*/ 

return (0) ; 

} 


FUNCTION  compl«x_math ( ) 

/*************************************************************** 
This  file  is  the  source  for  a  collections  of  subroutines 
written  by  Max  Blair  to  perform  simple  complex  operations. 
*******************************************«>********************/ 

♦include  <stdio.h> 

♦include  <math.h> 

♦define  ABS  (x)  (((x)<0)  ?  -  (x)  :  (x) ) 

♦define  EPS  (1.0e-20) 

♦define  BIG  (1.0e+20) 

/*  add  c**a+b  */ 

int  cadd(ar,ai,br,bi,cr,ci) 

float  arfai;/*  input  this  complex  number  */ 

float  br,bi;/*  input  this  complex  number  */ 

float  *cr, *ci; /*  return  this  complex  number  */ 

{ 

*cr  *»  ar+br; 

*ci  «■  ai+bi; 
return (0) ; 

} 

/*  add  c*a~b  */ 

int  caub(ar,ai,br,bi,cr,ci) 

float  ar,ai//*  input  this  complex  number  */ 

float  br,bi//*  input  this  complex  number  */ 

float  *cr,*ci//*  return  this  complex  number  */ 

{ 

*cr  «  ar-br; 

*ci  ■  ai-bi; 
return (0) ; 
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FUNCTION  compkxjoaathO 


} 

/*  multiply  c»a*b  */ 

int  cmult(ar,ai, br,bi,cr,ci) 

float  ar,ai;/*  input  this  complex  number  */ 

float  br,bi;/*  input  this  complex  number  */ 

float  *cr, *ci;/*  return  this  complex  number  */ 

{ 

*cr  *  (ar*br) - (ai*bi) ; 

*ci  *  (ar*bi) + (ai*br) / 
return ( 0 ) ; 

} 

/*  divide  c=a/b  */ 

int  cdiv(arf  ai,br,bi,cr,ci) 

float  ar,ai;/*  input  this  complex  number  */ 

float  br,bi;/*  input  this  complex  number  */ 

float  *cr,*ci;/*  return  this  complex  number  */ 

{ 

float  d/ 

d-  (br*br) +  (bi*bi) ; 
if (d<EPS) 

{ 

printf  (''An  division  by  complex  aero  in  cdiv\n") ,* 
exit  (0) ; 

) 

else 

{ 

*cr  *  ( (ar*br) + (ai*bi) ) /d; 

*ci  ■  ( (ai*br) - (ar*bi) )  /d / 

) 

return (0) ; 

) 

/*  transforms  complex  number  from  cartesian  to  polar  form  */ 
int  polar (ar,ai,bm, bp) 

float  ar,ai;  /*  input  cartesian  form  of  complex  number  a  */ 
float  *bm, *bp;  /*  return  polar  form  of  b,  mag  and  phase  (rad)*/ 
{ 

float  ftheta ( )  / 

*bm  -  sqrt ( (double) ( (ar*ar) +  (ai*ai) ) )  / 

*bp  -  ftheta (ar, ai) ; 
return  (0) ; 

) 

/*  transforms  complex  number  from  polar  to  cartesian  form  */ 
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FUNCTION  poweiQ 


int  cartesian  (am,  ap,br,bi) 

float  am, ap/  /*  input  polar  form  of  complex  number  a,  magnitude 
and  phase  (rad)  */ 

v  float  *br,*bi;  /*  return  cartesian  form  of  complex  number  b  */ 

{ 

*br  *  am*cos ( (double) ap) ; 

•  *bi  =  am* s in ( (double) ap) / 

return  (0) / 


/*  return  the  value  of  theta  (rad)  given  x  and  y  coordinates:  */ 
float  ftheta(x,y) 
float  x,y; 

{ 

float  pi, xtest, theta; 

pi=acos ( (double) (-1.0) ) ; 

xtest=f abs (  (double) (y*1.0e-05)  ); 

theta=pi/2 .0; 

if (y<0. 0) theta® (-pi) /2. 0; 

if  (fabs  ( (double)  x)  Oxtest)  return  (theta)  ; 

theta*atan  ( (double)  (y/x) ) ; 

if (x<0 . 0) theta»theta+pi; 

return (theta) / 

} 

/*  absolute  value  of  complex  number  in  cartesian  form  */ 
float  cmag(ar,ai) 

float  ar,ai;  /*  input  cartesian  form  of  complex  number  a  */ 

{ 

float  mag;  /*  return  polar  form  of  b,  magnitude  and  phase  (rad)*/ 
if  (  ABS (ar) >BIG  (|  ABS (ai) >BIG  ) 

{ 

printf ("potential  error  in  cmag  ar:  %e  ai:  %e\n#, ar, ai) ; 
if(  ABS (ar) >BIG  )mag  ■  ar; 
if (  ABS(ai)>BIG  )mag  *  ai; 
return (mag) ; 

) 

if (  ABS (ar) <EPS  &&  ABS (ai) <£FS  ) return (0.0) ; 
mag  ■  sqrt ( (double) ( (ar*ar) + (ai*ai) ) ) ; 
return (mag) ; 

) 
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FUNCTION  qusdrilateralO 


float  power  (x,i) 

/*  carry  out  the  operation  xAi  where  x  is  real  and  i  is  integer  */ 
int  i? 
float  x/ 

{ 

int  j; 
float  p / 

if ( i“0 ) return  ( 1 ) /  « 

if  (i=l)  return  (x) ; 

P*x;  t 

for (j»2/ j<«i; j++) 

{ 

p«p*x; 

} 

return (p) ; 

} 


FUNCTION  quadrilateral ( ) 

/*  input  the  coordinates  in  a  counter-clockwise  order  */ 

♦include  Ml.  define" 

♦include  <math.h> 

♦include  <stdio.h> 

♦include  *dl.  structure*' 

int  quadrilateral ( odat , xl , y 1 , x2 , y2 , x3 , y 3 , x4 , y 4 , area , xc , y c ) 

FILE  *odat; 

float  xl,x2,x3,x4; 

float  yl,y2,y3,y4; 

float  *area,  *xc,  *yc; 

{ 

float  al,a2,a3,a4; 
float  bl,b2,b3,b4; 
float  01,02/03; 

al  «  (-xl+x2+x3-x4) /4.0; 
a2  «  (  xl-x2+x3-x4) /4. 0; 
a3  *  (-xl-x2+x3+x4) /4.0; 
a4  ■  (  xl-x2+x3-x4) /4. 0; 
bl  »  (-yl-y2+y3+y4) /4.0; 
b2  -  (  yl-y2+y3-y4) /4.0; 
b3  m  (-yl+y2+y3-y4) /4. 0; 
b4  -  (  yl-y2+y3-y4) /4.0; 
cl  -  (  al*bl-a3*b3) ; 
c2  «  (  al*b2-a2*b3) ; 
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c3  <*  (  a2*bl-a3*b2) ; 

*area  =  4.0*cl/ 
if (*area<=0.0) return (1) ; 

al  =  {  xl+x2+x3+x4) /4 . 0/ 
a2  <=  (-xl+x2+x3-x4) /4.0; 
a3  =  (-xl-x2+x3-fx4) /4 . 0; 
a4  *  (  xl-x2+x3-x4) /4.0; 
bl  =  (  yl+y2+y3+y4) /4.0; 
b2  =  (~yl+y2+y3-y4) /4 . 0; 
b3  *  (-yl-y2+y3+y4) /4.0; 
b4  *  (  yl-y2+y3-y4) /4.0; 

*xc  =  (  4.0* (al*cl)  +  4.0* (a2*c2)/3.0  +  8.0*(a3*c3)  )/{*area) 
*yc  «  (  4.0* (bl*cl)  +  4.0*  (b2*c2)/3.0  +  8.0*(b3*c3)  ) / (*ar©a) 

return (0) / 

} 


FUNCTION  inJL.in«M 

♦include  <stdio.h> 

int  in-line (idat, line) 

char  line [ ] / 

FILE  *idat; 

{ 

char  cl; 
int  j; 

cl»getc(idat) ; 

if  (cl— EOF) 

{ 

line [0]“' \0' ; 
return (1) ; 

} 

elae  if  (cl— 10) 

{ 

line[0]-'\0'; 
return  (0) ; 

) 

else 

{ 

line  [0) -cl; 
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end  of  doOUet  Isttice  source  code 


} 

while (  (line! j]*getc (idat) )  !-  10  ) 

{ 

if  (line  [j]  “EOF) 

{ 

printf ("\nEND  OF  FILE  ENCOUNTERED  IN  in-line  () \n") ; 
exit  (0) / 

> 

j++; 

} 

line[  j]«'  \0'  / 
return (0) ; 

} 


APPENDIX  B 


The  Doublet  Lattice  Program  Input  File 


digram 


BLAIRCRAFT  2100  ATTACK  FIGHTER 


6.0 

0.5 

1.00 

s 

0.0  0.0 
12.0  0.0 
12.0  12.0 

0.0  12.0 
3 
3 


characteristic  length  (b) 

Mach 

reduced  frequency  wb/U 

s:  symmetric  a:  anti -symmetric  n:  no  symmetry 
x  and  y  coord  of  inboard  leading  edge 

x  and  y  coord  of  inboard  trailing  edge 

x  and  y  coord  of  outboard  trailing  edge 

x  and  y  coord  of  outboard  leading  edge 

number  of  chordwise  cuts  (discretized  x) 
number  of  spanwise  cuts  (discretized  y) 


COMMENTS : 


Pdimensional  pressure 
Cpnon-dimensional  pressure 
rhoair  density 

P  -  rho*UA2*Cp 

For  typical  "non-dimensional"  input,  set  U*b»l. 
The  output  is  interpreted  accordingly. 

Only  Cp  is  printed  out  for  each  box. 


bc_INPUT 

Boundary  condition  (upwaah)  input  for  doublet  lattice: 

flag  constant  x  power  y  power 

1  -1.0  0  0 

0  -1.0  1  0 

0  -1.0  0  1 
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2 

1 

0 


0 

1 

2 


0  -1.0 

0  -1.0 

0  -1.0 

end  of  data 

interpretation : 
w(x,y)  ■  aOO  +  al0*x  +  a01*y  +  a20*x*2  +  all*x*y  +  ai)2*yA£ 

instructions : 

Only  data  with  a  "1"  in  the  first  coluxnh  viil  be  cohsiddfed 
data. 

Replace  the  *1"  with  a  "0"  to  ignore  any  data. 

A  line  which  begins  with  an  "e*  will  terminate  the  input ; 
There  must  be  at  least  one  line  which  bfegin M  with  an  *e". 


APPENDIX  C 


The  Doublet  Lattice  Program  Output  Listing 


Auxiliary  runtime  data  placed  in  file  [dl. TRASH] 

Aerodynamic  influence  coefficients  placed  in  file  [dl.AIC] 
MAXBOX:  400  MAXDIM:  160Q00 
Begin  input 

Input  data  will  be  read  from  file  [dl. INPUT] 

TITLE: 

text}  [BLAIRCRAFT  2100  ATTACK  FIGHTER] 
characteristic  length:  €.000000 
Mach:  0.500000 
reduced  frequency:  1.000000 
Assume  symmetry  about  the  x  axis 

text:  {  0.0  O.Ox  and  y  coord  of  inboard  leading  edge] 
inboard  leading  edge:  x  0.000000  y  0.000000 
inboard  trailing  edge:  x  12.000000  y  0.000000 
outboard  trailing  edge:  x  12.000000  y  12.000000 
outboard  leading  edge:  x  0.000000  y  12.000000 
number  of  boxes  in  the  x  direction:  3 
number  of  boxes  in  the  y  direction:  3 
Input  complete 

The  wing  area  used  to  non-dlmensionalize  lift  is:  4.0000e+00 
The  wing  centroid  is  at  x:  1.000000  and  y:  1.000000 
Begin  discretizing  the  wing. 

Discretization  is  now  complete, 
reduced  freq  (k) :  1.0000e+00 
beta*2  <1-MA2) :  7.5000e-01 

0123456780 

12345678 
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Boundary  condition  data  will  be  read  from  file  [dl. INPUT] 
reduced  frequency  (k) :  1.0000e+00 
[1]  a:  -l.COOOOO  px:  0  py:  0 

compute  the  upwash  at  9  control  points 


Upwash 
0  Real 

1  Real 

2  Real 

3  Real 

4  Real 

5  Real 

6  Real 

7  Real 

8  Real 


specified  and  pressure 


[w] 

[w] 

[w] 

[Wj 

[w] 

[w] 

CW] 

[w] 

[w] 


0.0000e+00  Imag[w] : 
O.OGOOe+QO  Imag[w] : 
0.0000e+00  Imag[w]: 
0.0000e+00  Imagtw] : 
0 . 0000e+00  Imag [w] : 
0.0000e+00  Xmag[w]  : 
O.OOQOe+OO  Imag[w]: 
0.0000e+00  Imag[w]: 
0.0000e+00  Imag[w]: 


will  be  confuted 
-1. 0QQ0e+00 
-1. 0000e+00 
-1. 0000e+00 
-1 . 0000e+00 
-1. OQOOe+QO 
-1. 0000e+00 
-1. 0000e+00 
-1 . 0000e+00 
-1. 0000e+00 


Solve  the  complex  problem  {w}  *  [D] {p> 


Cp  PRESSURE  COEFFICIENTS  (P-0 . 5*rho*U*2*Cp) 
box  #  (real  Cp)  (imag  Cp)  (box  area) 


0  -5. 4900e-01 

1  -3.88626+00 

2  '-3. 8736e+00 

3  ~5.9144e“01 

4  -3. 6405e+00 

5  -3. 6234e+00 

6  -5 . 8286e-0l 

7  >2.8983e+00 

8  -2 . 8893e+00 


6.2682e+00  1.6000e+01 
2 . 4495e+00  1.6000e+01 
1 . 1745e+00  1 . 6000e+01 
5. 8092e+00  1.6000e+01 
2 . 1530e+00  1.6000e+01 
1.0281e+00  1 . 6000e+01 
4.5474e+00  1.6000©+01 
1 . 4663e+00  1.6000©+01 
7 . 1186e-01  1 . 6000e+01 


The  characteristic  length  b:  The  non-dimensional  ving  area  used 
to  non-dimensionalise  lift  is:  4.0000e+00 
Lift  Coefficient  «  C_L*q*A 

THE  COMPLEX  WING  LIFT  COEFFICIENT  (C_L)  IS: 

[  (-2.50386+00)  +  (  2. 8433e+00) i] 

MAGNITUDE:  (  3.7901e+00)  PHASE:  (  131.3471)deg 


end  of  data 
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U.l.  CoMroMet  Mating  Ofttc*  MI-11? 


